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LIFTING-LINE THEORY FOR A WING IN NON-UNIFORM FLOW* 


BY 
THEODORE VON KARMAN anp HSUE-SHEN TSIEN 
California Institute of Technology 


1. Introduction. Prandtl’s theory of the lifting line gave the answer to most of the 
questions in the aerodynamic design of airplane wings. Thus the three-dimensional 
wing theory became a standard tool of airplane designers. One restriction involved 
in the conventional wing theory is the uniformity of the undisturbed flow in which 
the wing is placed. Now there are many important cases which do not satisfy this 
condition. For instance, in the case of a wing spanning an open jet wind tunnel, the 
velocity of the air stream has a maximum at the center of the jet and drops to zero 
outside of the jet. Another example is the problem of the influence of the propeller 
slip-stream on the characteristics of the wing. Here the higher velocity of the propeller 
slip-stream makes the application of the Prandtl wing theory difficult. Such cases led 
several authors to investigate the problem of a wing in non-uniform flow. Some in- 
vestigators found a satisfactory solution of the problem for the case of “stepwise” 
velocity distribution. In this case the flow in regions of uniform velocity can be deter- 
mined by using Prandtl’s concepts with additional continuity conditions at the bound- 
aries between such regions. On the other hand, the problem of a continuously varying 
velocity field seems to need an appropriate treatment. K. Bausch! has tried to modify 
the Prandtl theory for the case of small inhomogeneity in the air stream; however, 
besides the restriction of slight deviation from uniform flow, his method encounters a 
further difficulty in estimating the error introduced by the approximations. The 
seriousness of this difficulty becomes evident when one tries to compare the results of 
Bausch with that of F. Vandrey.? Vandrey considers the problem with variable 
velocity as the limiting case of a wing in a stepwise velocity field, and his result seems 
to differ from that of Bausch. Recently R. P. ISaacs* has investigated the same prob- 
lem, but the authors have not yet had the opportunity to study his work. 

It seems to the authors that a general and more satisfactory solution for the flow 
of a wing in a non-uniform stream can be obtained by studying the three-dimensional 
problem anew in this generalized case, introducing the modifications of Prandtl’s 
fundamental concepts. The first fundamental concept is the following: the span of 


* Received September 27, 1944. 
1K. Bausch, Auftriebsverteilung und daraus abgeleitete Grissen fiir Tragfliigel in schwach inhomogenen 


Strémungen, Luftfahrtforschung, 16, 129-134 (1939). 
2F, Vandrey, Beitrag sur Theorie des Tragfliigels in schwach inhomogener Parallelstrimung, Zeit- 


schrift f. angew. Math. u. Mech. 20, 148-152 (1940). 
3 R. P. Isaacs, Airfoil theory for flows of variable velocity, abstract in Bulletin of the American 


Mathematical Society, 50, 186 (1944). 
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the wing is sufficiently large compared with the chord so that the variation of the 
velocities in the spanwise direction is small when compared with the variation of the 
velocities in a plane normal to the span; then the flow at each sectional plane per- 
pendicular to the span can be considered as a two-dimensional flow around an airfoil. 
The only additional feature for the flow in this sectional plane is the modification of 
the geometrical angle of attack, as defined by the undisturbed flow, on account of the 
so-called induced velocity. The second fundamental concept of Prandtl is the replace- 
ment of the wing by a lifting line having the same distribution of lifting forces along 
the span as the wing. This concept, 

z with the additional assumption that 
i, w the disturbance caused by the lifting 
t av line is small, i.e., that the wing is 
UY+u lightly loaded, makes the calculation 


of theinduced velocity relatively sim- 


po) YW ple. In this paper the authors will 

ei study the flow around a lightly loaded 

—PF Uy, Z) ha sib biriatead lifting line placed in a parallel stream 
Ml a" whose velocity is perpendicular to the 





span (Fig. 1) and is assumed to vary 
in both directions normal to the flow. 
Due to the rather complicated char- 
acter of the flow, the usual concept of 
the picturesque system of trailing 
vortices encountered in Prandtl’s 
wing theory is not very useful here. 
A method, which is mathematically 
Fic. 1. Lifting line in a non-uniform flow. more convenient, has to be adopted. 
This method has already been used 
by the senior author‘ in explaining the similarity between Prandtl’s wing theory and 
the theory of planning surfaces. After the general theory is formulated, the problem 
of minimum induced drag will be considered. Finally a general expression for calculat- 
ing the induced drag of a wing in a stream of varying velocity will be presented. 

Of course, the complete solution of the problem of a wing in a non-uniform stream 
requires a knowledge of the “section characteristic” or the two-dimensional properties 
of the airfoil sections of the wing. If the velocity of the main stream is varying only 
in the direction of the span, the required section characteristics are those of an airfoil 
in a two-dimensional uniform flow, and are common knowledge in applied aerodynam- 
ics. However, if the velocity of the main stream is also varying in a direction perpen- 
dicular to the span and to the velocity itself, the required section characteristics are 
those of an airfoil in a two-dimensional non-uniform flow. Such flow problems have not 





yet been studied extensively.® 
2. General theory of a lifting line. Let the x-axis be parallel to the direction of 
the main flow, the y-axis coincide with the lifting line and the z-axis be normal to the 





4 Th. von K4érm4n, Neue Darstellung der Tragfliigeltheorie, Zeitschrift f. angew. Math. u. Mech. 15, 


56-61 (1935). 
5 H.S. Tsien, Symmetrical Joukowsky airfoils in shear flow, Quart. Appl. Math. 1, 130-148 (1943). 
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lifting line (Fig. 1). If p is the pressure, p the density, and 2, v2, vs the components of 
the velocity, the dynamical equations for the steady motion of an inviscid, incompres- 
sible fluid without external forces are 





Ov Ov; Ov, 1 dp 

V9 Oy Qe ee es (1) 
Ox oy 02 p Ox 
Ove Ove OV 1 Op 

_. ea + ara oa V3  - <= om, (2) 
Ox oy 02 p oy 
Os OV3 03 1 @d 

Pops a ee Pe ns. a (3) 
Ox oy p oz 


The equation of continuity is 


Ov Ove Ov3 


— = (), 4 
Ox Oy Oz (4) 


Equations (1) to (4) constitute a system of four simultaneous equations for the four 
unknowns 2), %, ¥3 and p. 

For the particular problem of a lightly loaded lifting line, the velocity components 
can be expressed in the following forms: 


1=U+4u, (5); Ve = 2, (6); v3 = W. (7) 


Here u, v, w are the velocity components due to the presence of the lifting line and U 
is the main stream velocity assumed to be a function of y and z but independent of x. 
Since the lifting line is assumed to be lightly loaded, u, v and w are small compared 
with the main velocity U. By substituting Eqs. (5) to (7) into the dynamical equa- 
tions and neglecting higher order terms, a set of linear equations for u, v and w is 
obtained. Thus 

aU aU aU 1 ap 








U—+r7—+w aS Sy (8) 
Ox Oy Oz p Ox 
Ov 1 0 Ow 1 0 
U—=-— =, (9); U =z=—— “ f (10) 
Ox p oy Ox p Oz 


Then the equation of continuity becomes 


Ou Ov Ow 


—+—+ 
Ox oy Oz 


I 
= 





(11) 


If Eqs. (8), (9) and (10) are differentiated with respect to x, y and z respectively and 
the results added, the sum can be simplified by using Eq. (11) and can, finally, be 


1 ap 0/1 op 0/1 op 
= t+ a(— 7 )+2(— 2) -0. (12) 
U? dx? dy\U? dy 0z\ U2? dz 


This is now an equation for the pressure p only and can be used conveniently as the 
starting point of the solution, If the pressure of the undisturbed main flow is chosen 


written in the form 
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as the reference pressure and set equal to zero, one of the boundary conditions to be 
satisfied by p is 
| | 
p=0, for |x| @, | y|— @, or || —> @, (13) 
The condition at the lifting line, or y-axis, is that the lifting force is represented by a 
suction force on the “upper surface” of the lifting line and a pressure force of equal 
magnitude on the “lower surface” (Fig. 2). Hence the pressure p must satisfy the 


following expressions 





p | pdx = 3l(y), for z= oh 0, (14) 
- 
le and 
>< 
5) 7 a : 
| x= sl(y f = 0, 1 
PRESSURE ON LOWER y J mas 5l(y), or 2 (15) 
SURFACE Lhe 


where /(y) is the lift per unit length of the 





t x lifting line at the point y. Furthermore, on 
account of the symmetry of the flow, 
Wve 
SUCTION OW UPPER : p=0 for z=0, |x| >. (16) 
SURFACE =I ) x ‘ 
im To solve Eq. (12) together with the 
Ile boundary conditions given by Eqs. (13) to 


(16), the Fourier integral theorem can be 
used to build up the solution of the problem 
from the elementary solutions of Eq. (12) 
Fic. 2. Representation of lift as pressure forces of the form 

acting on the two “surfaces” of the lifting line. P(y, 2, ) cos dx. 


The equation to be satisfied by P is 


a/1 aP a/1 oP 
r oe 1h et Se lee’, 17 
;' sla 5 )t =a =) (17) 


To determine P uniquely, it is convenient to impose the following conditions 


P=0, for |vj/ 7, |z|}> 2, (18) 
P=-—3(y) for z= +0, (19) 
P=il(y) for z= —0. (20) 


The required solution for p can then be written as 


1 ee) 
p= —f cos AxP(y, z, A)dd. (21) 
0 


T 


By substituting Eq. (21) into Eqs. (9) and (10), the “induced velocities” v and w 


are obtained; 


,..4 * sindkx @ 
v(x, y, z) = v(0, y, 2) — f — — P(y, 2, d)dd, (22) 
pU rvJo d Oy 
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sindkx 0 
— P(y, z, )dd. (23) 
r Oz 





«<8 o 
w(x, y, 2) = w(0, y,z) — — —f 
pU T/0 


Because the integrals are odd functions of x, the following relations hold for velocities 
far ahead of the lifting line and far behind the lifting line: 


+ [o(— ©, Y, z) + v( 0, y; z) | _ v(0, y; z), 4[w(— o,Y, z) + w(o, y; z) | ” w(0, ys z). 


However, it is evident that the induced velocities far ahead of the lifting lines must 


be zero. Hence 
v(0, y, 2) = 3n(0, y, 2), (24); w(0, y, s) = }w(, y, 2). (25) 


The induced velocities v and w at the lifting line are then one-half of those far down- 
stream. This is in accordance with the usual wing theory based upon the concept of 
trailing vortices. 

One meets an apparent difficulty if the x component of the induced velocity is cal- 
culated; integration of Eq. (8) with respect to x furnishes the x-component of the 
induced velocity: 


1 1 oU f? 1 ov f* 
om ——p—- — — vdx — — — wdx. (26) 


pl U dyJ_. S 88 VP ire 


Since p tends to zero, v and w tend to finite quantities as x tends to infinity, and u in- 
creases indefinitely as x tends to infinity. This is in contradiction to the assumption 
of small disturbances introduced at the beginning of the present investigation. How- 
ever, it is believed that this difficulty does not prevent the application of the theory to 
practical cases, since the apparent large value of the u component is due to the dis- 
tortion of the variable main stream by the induced cross flow and the infinite value 
for x—> is due to the linearization of the differential equations. Some further re- 
marks on this point are given in Section 4. 

3. Conditions far downstream. For the application of the lifting-line theory to the 
wing problem, the quantity of primary interest is the z component of the induced 
velocity at the lifting line. The simple relations given by Eqs. (24) and (25) suggest 
a possible simplification of the calculation by considering conditions far downstream, 
or the “Trefftz plane” according to the terminology of the conventional wing theory. 
To abbreviate the notation, we let 


vo = v(0, y, 2), Wo 


w(0, y, 2), \ (27) 


v1 = (0, y, 2), wi = w(~, y, 2). 


Then, according to (24) and (25), 10 = 401, wo = 4m. Therefore, Eqs. (22) and (23) give 








eer * sindx @ 

v1, = — — lim —f — P(y, z, \)dx, 
pU z~e2 0 r Oy 
:. 2 f*ck @ 

w; = — — lim —f — P(y, z, Add. 
pU zo 0 PN Oz 


Let us consider P(y, 2, \) as a regular function of A; then 
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oP 
P(y, 2, +) = P(y, 2, 0) + | + 
On 


A=0 


By using the variable =x, the expressions for v; and w can be rewritten, 


1 2 f* sint @ t (oP 
1, = — — lim —f — | Po. s, 0) + — (—) atc Ja 
pU z308 Tr J 9 t oy x \ OX /y=0 


ie * sint @ . t (oP 
4 2 = = lim — spores | Po. Z, 0) + ~(=) te: | dt. 
0 x Or /r=0 


pU z30 TJ $ OB 








At the limit, only the first terms of the integrands are significant, and furthermore 
2 f* sint’ 
_— — dt = 1, 
T 0 t 
1 


0 i @ 
n= ———P(y,2z,0), (28); w, = —— —P(y,2z,0). _—.(29) 
pU dy pU dz 


Hence 
1 


Equations (28) and (29) simplify the problem of calculating the induced velocities 
at the Trefftz plane considerably. In fact, by introducing a “potential function” ¢ de- 
fined by the relation 


o(y, 2) = — Ply, z, 0), (30) 


the problem can be formulated as follows: the differential equation to be satisfied 
by ¢ can be deduced from Eq. (17) by setting \=0; thus 


afi a/1a 
<( “+ —( =) =0. (31) 
ay\U? ay) dz\U? az 


The boundary conditions to be satisfied by ¢ are 








@=0 for | y|-> 2, |z|— @, (32) 
@=l(y)/2 for z=+0, (33) 
@=—IU(y)/2 for z= —0. (34) 
Then 
1 d¢ 1 d¢ 
je —) (35); {ee —* (36) 
pU day pU az 


By substituting Eqs. (35) and (36) into Eq. (31), one has 


re) V1 te) WwW 
—(*) + ~(=) = 0. (37) 
dy \U dz \U 


This equation has a very simple physical meaning. Since 1; and w; are considered to 
be small quantities, the ratios v1/U and w:/U are the angles of inclination, 6 and y, 
of the stream lines with respect to the zx and xy planes. Consider parallel planes 
perpendicular to the x-axis and dx apart (Fig. 3). If the width of the stream tube 
at the section x is 6,, then at the section x+dx, the width of the stream tube is 
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5,[1+dx 08/dy]. If the height of the stream tube at the section x is 6,, then at the 


section x +dx, the height of the stream tube is 6,[1-+dx dy/0z]. The total increase in 
the cross-sectional area of the stream tube from x to x+dx is then approximately 


x 
5a 
x a |h , 
- 7+ 26, “a 8, (14 Ex) 
, 
Seer oes 
+7 4 213- 


are 


























Fic. 3. Stream tube far downstream from the lifting line. 


Oop ay 
5): ( + ~) dx. 
Oy 02 


Now at the Trefftz plane, the flow field can be considered as settled into a uniform 
condition; i.e., the pressure is constant in the x-direction. Hence, the velocity of the 
flow along any stream tube is constant. Then the cross-sectional area of the stream 
tube must be also constant. Therefore, 


op oY 

re + a 0, 

dy O02 
which is simply Eq. (37). From this point of view, Eq. (37) is really the equation of 


continuity, simplified under the conditions prevailing at the Trefftz plane. 
On the other hand, ¢ can be eliminated from Eqs. (35) and (36). The result is 


* (Un) — —(Um) = 0 (38) 
a v =— eo w = as x“ 
as Oy ; 


This equation can be considered as the modified vorticity equation. It actually holds 
for all values of x under the approximation assumed in the present investigation. This 
can be seen in the following way: since U is a function of y and z but independent of x, 
Eqs. (9) and (10) can be written in the form 

a 1 ap o. 1 ap 

—Uvse=z -—- — —» —Uw= -— — - 

Ox p oy Ox p oz 
By differentiating the first equation with respect to z and the second equation with 
respect to y and then subtracting, the result is 
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—_ U U = 0 
x5 wo - 5 um] =o. 


0 
(Uv) — = (Uw) = a function of y and z 
. 


Thus 
0 


But for points far upstream, or for x = — ©, v and w vanish; therefore the function 
of y and z on the right of above equation must be identically zero. Hence for all 


values of x, 


— it~ fw w) = 0. (39) 
0z oy . 
It should be noted here that Eqs. (37), (38) and (39) are obtained without any 
reference to the lifting line and hence they are true for more general cases. However, 
the complete determination of v; and w; requires a knowledge of the relation between 
the induced velocities and the lift on the wing. This relation depends upon the type 
of lift distribution. For the particular case of a lifting line, this relation is supplied 
by Eqs. (33) and (34) 
Equation (37) can be identically satisfied by introducing the “stream function” 
WY defined by 
71= U Z “14=> => om: (40) 
Oz oy 


Then Eq. (38) gives the differential equation for y: 


) Oy fe) oy 
— (us ~) + —(u ~) _" (41) 
Oy Oy Oz 0z 


Both Eq. (31) and Eq. (41) reduce to the Laplace equation for the conventional wing 


theory when U is a constant. 

4. Minimum induced drag. The induced downwash angle at the lifting line is 
equal to wo/U or 4w,/U, according to Eq. (25). Therefore, Eq. (36) gives the down- 
wash angle at the lifting line as [1/2pU*](0@/dz),~0, and the induced drag D; can then 


be expressed as 


0g 
on = - =f lo(y, + 0) — o() , — 0)] —(° =) iy- > f ~~ ds. (42) 
02/ 2=0 Cc 2 02 


The first integral is evaluated across the span of the lifting line. The second integral 
is calculated along a contour following the upper and lower “surface” of the horizontal 
strip shown in Fig. 4. Since ¢—0 for points far from the lifting line, the contour in- 
tegral can be transformed into an area integral by Green’s theorem, and 


1 d¢ 1 d¢ 
- + — (= )} dyd. (43 
-=-ff ee 72 o) 6z\U ar — 


This integral extends throughout the region outside of the lifting line. Since ¢ satisfies 
the differential equation (31), Eq. (43) reduces to 
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»- SiGe Gah wo 


Therefore, the induced drag is represented by the kinetic energy corresponding to the 
velocity components 2 and w; at the Trefftz plane. It is seen that the u component 
of the velocity does not appear in the expression for the induced drag. This is due to 
the fact that the increase of u with increasing x does not represent a real acceleration 











z A 
UW, 
toy 
PROJECTION OF LIFTING LINE 
ON TREFFTZ PLANE 
1, = 2 —=> —>y 





Fic. 4. Contour integration in the Trefftz plane. 


of a fluid element in the x direction. Rather, it is due to the fact that the cross flow 
transports fluid elements from regions of lower main velocity to regions of higher 
main velocity and vice versa. This is in accordance with the modified continuity 
equation (37) which clearly indicates that the cross section of the individual stream 
tubes has a definite limiting value for x, and therefore the velocity component 
in the direction of the stream tube tends to a finite value. 

The problem of minimum induced drag requires the determination of the mini- 
mum of D; as given by Eq. (44) together with the condition that the total lift LZ 


remains fixed. Thus 
L= fitay = f [o(y, + 0) — o(y, — 0) |dy = — f gods = constant. (45) 
o 


By using the method of Lagrange’s multiplier, the above problem can be reduced to 
that of finding the minimum of D;+K/pL, where K is a constant. Hence, 


K 
6D; + —aL = 0. (46) 
p 


The variation of the induced drag can be obtained from Eq. (44), 


1 1 dg 1 06 1 d¢ 1 0 
=f fe pas 
p U dy U day U dz U Oz 


However, ¢ must satisfy the differential equation (31) ; thus 


6D ~ff (—° v ) + : —( ps ie) lydz 4 %, d. 
i, — yaz = uname 
ts U? oy ° U? dz fa Cc U? dz wits 
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éL = — f seas. 
e 


By substituting these results into Eq. (46), the condition of minimum induced drag 
is obtained in the form 

1 | 

pJc\U? az 


The variation of 6¢ on the lifting line is arbitrary; therefore the minimum induced 
drag is given by the condition that the induced downwash angle must be constant 
along the span. If the main stream velocity U is constant, the above condition is 
reduced to the requirement of constant downwash. This is in agreement with the 
well-known result of Prandtl’s wing theory. 

5. Flow with velocity varying in the direction of span only. If the stream velocity 
varies only in the y direction, i.e., in the direction of the wing span, the calculation 
of induced velocity and induced drag can be simplified with the aid of characteristic 
functions connected with the differential equation for the potential function ¢. In this 


10 


On the other hand, 


case Eq. (31) becomes 








dU 
ap a dy 0d 
- — 2?—— — = 0. (48) 
dy? az? U_ ody 


To satisfy the boundary condition given by Eq. (32), @ is expressed by the following 





integral 
o(y, 2) -{ f(r\)e™*V(y) dr (49) 
0 
for z>0. f(A) is an unknown function to be determined. For 2<0, 
: o(y, 2) = — o(y, — 2). (50) 
By substituting Eq. (49) into Eq. (48), the differential equation for Y,(y) is ob- 
tained, ; 
dU 
a°*V, dy dV, 
— 2—— — +), = 0. ($1) 
dy? L dy 


This equation will determine Y,(y) uniquely if proper normalizing and boundary con- 


ditions are imposed. 
At the span, the condition (33) must be satisfied. Thus 


Uy) anni ; 

—— =| f)Fi(y)dd. (52) 
2 0 

This relation can be considered as the equation for determining f(A) with the given 


lift distribution /(y). For example, in the case of constant stream velocity U or 
Prandtl’s case, Y,(y) is a trigonometric function and therefore f(A) can be deter- 
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mined easily by means of Fourier’s inversion theorem. Equation (50) shows that with 
f(A) so determined, the condition (34) will be automatically satisfied. 

The downwash velocity wo at the wing can then be easily calculated by using 


Eqs. (25), (36) and (49). The result is 


1 i) 
wo(y, 0) = — at A(A)Va(y) dr. (53) 
pu #90 


“ o(y, 0 
a J I(y) ers 
rn U 


Therefore, in terms of Y,(y), the following general expression for the induced drag 


The induced drag D, is given by 


is obtained: 


D; = i ——ay f F(NVa(y)dd é nf(n)¥ (y)dn. (54) 


» pl 


Thus the problem of calculating the induced drag with a given distribution of lift 
I(y) is reduced to the problem of solving the integral equation (52) for f(A) and then 
evaluating the integral given by Eq. (54). 

If the chord c, the geometrical angle of attack a and the slope k of the lift coeffi- 
cient are given instead of the lift distribution /(y), then 


wo(y, 0) 
(y) = 3pU chap eh (55) 


Thus Eq. (52) is replaced by the following equation 


1 2) is) 
tpU*ck {a _ —— f Or Fnin)an} = f TNVi(y) aA, 
2pU? Jo 0 


° k 
tpU*cka = f (1 + <n) f0)¥a(y)da. (56) 


This is now the integral equation for f(A). When f(A) is determined, the induced drag 
D; can be again calculated by using Eq. (54). 


or 








ON THEODORSEN’S METHOD OF CONFORMAL MAPPING 
OF NEARLY CIRCULAR REGIONS* 


BY 
S. E. WARSCHAWSKI 


Washington University 


1. Introduction. In determining the complex velocity potential of the two-dimen- 
sional flow around an airfoil, one is lead to the problem of finding the analytic func- 
tion which maps the exterior of a circle conformally onto that of a “nearly circular” 
contour. T. Theodorsen developed a method for the practical computation of this 
mapping function, a method which was later elaborated on in a joint paper by 
Theodorsen and I. E. Garrick.! Theodorsen reduces the problem of determining the 
mapping function to the solution of a certain non-linear integral equation which 
then is solved by successive approximations. In both papers examples of wing sections 
of airplanes are calculated demonstrating the use of the process and the rapidity with 
which it converges. However, the validity of the method from a mathematical point 
of view, such as the proof of the convergence of the successive approximations, is not 
discussed. The present paper is an attempt to supply such a discussion. Simple con- 
ditions on the nearly circular contour (essentially involving the tangent angle and the 
curvature) are established which insure the convergence of the process. The absolute 
value of the difference between the mapping function and the successive approxima- 
tions is estimated. These estimates serve both to prove the convergence and to ap- 
praise the accuracy of the approximation. The analogous problem for the derivative 
of the mapping function is treated. (The derivative of the mapping function enters 
in the computation of the velocity and pressure distribution on the surface of the 
wing.) Finally, conditions are discussed under which the map of the circle by means 
of the successive approximations is star-shaped. 

Although Theodorsen’s method is of particular importance in the theory of air- 
foils, it represents the solution of a general problem in conformal mapping. For this 
reason all results of the present paper are derived for the “standard” case where the 
interior of a circle about the origin is mapped onto the interior of-the nearly circular 
contour containing the origin under preservation of the positive line element at the 
origin. However, all results obtained remain the same for the mapping function of 
the exteriors and for a different normalization of the mapping function (see §3). 

Sections 2—8 contain the actual results and proofs of the paper. To simplify the 
presentation some auxiliary results used in the text are listed in §9. 

2. Theodorsen’s integral equation and the successive approximations. Let C bea 
simple closed curve represented in polar co-ordinates by the equation p=p(@) 
(0027), where p(@) is absolutely continuous? and for some e (0<€<1), 


* Received May 9,°1944; presented to the American Mathematical Society, November 26, 1943. 
1T. Theodorsen, Theory of wing sections of arbitrary shape, NACA Tech. Rep. No. 411 (1931); 
T. Theodorsen and I. E. Garrick, General potential theory of arbitrary wing sections, NACA Tech. Rep. 


No. 452 (1934). 
2 A function g(@) is absolutely continuous in an interval if its derivative g’(6) exists for all @ of this 
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a 
— Ss p(9) S al , 2.1 
pq 3M) Sal + 6) (2.1) 
a being a positive constant, and 
p’ (8) 
Se. (2.2) 
p(8) 





Any curve C satisfying these conditions will be called a nearly circular contour. 
Let us suppose that the function w=f(z) maps the circle | z| <1 conformally onto 
the interior of C, and that f(0) =0, f’(0) >0. The function 


fe) _ | f@) 


z 


+i arg) 
z 











F(z) = log log , (2.3) 








which is defined as the real-valued log f’(0) when z=0, is single-valued and analytic 
for | z| <1 and continuous for | s| <1. For z=e* we write arg [f(e*)e—*#] =0(¢) —¢, 
and therefore 


F(e'*) = log p[0(¢)] + i(0(¢) — ¢). (2.4) 


Hence 
1 . t 
0¢) —¢ = — —f {log p[6(@ + #)] — log p[0(@ — #)]} cot su (2.5) 
T/ 0 


(The term arg [f(z)/z]..o=arg f’(0), which should be added to the integral on the 
right, is zero.) Thus the function F(e*) and hence f(z) may be found by solving this 
integral equation for 0(@). The existence of a continuous solution of this integral equa- 
tion is assured by Riemann’s mapping theorem. This solution is also umique as is 
shown in §9(a). In order to compute the solution we follow Theodorsen and form the 
successive approximations 


8o(o) = ?, 


1 * t 
6.(¢) -¢ = — —f {log p[0,-1(@ + #)] — log p[0,-1(@ — 2)]} cot — dt (2.6) 


(n = 1,2,-+-). 


The functions 0,(@) are continuous for 0 So S27; in fact, they are absolutely continuous 
and the squares of their first derivatives are integrable (for the proof of this, see 
§9(b)); 6,(¢) —@ is a conjugate function of log p[8,-1(¢) |. 

We shall show that the sequence 0,(@) converges uniformly to 0(@) as n—«. Hence, 
also log p[@,(¢)] converges uniformly to log p[@(¢)] as n>, so that the functions 


Fy(e*) = loga, F,(e**) = log p[n1(¢)] + i(6.(@) —¢) (m2 1), (2.7) 


may be used to compute F(e**) with any desired degree of accuracy. 
Let F,(z) denote the function which is analytic for || <1 and assumes the bound- 
ary values F,(e**) for | z| =1. By the principle of the maximum modulus the uniform 





interval except possibly for a set of Lebesgue measure zero and if JS2z'(0)d0 = g(b) —g(a) for every a and b 
of this interval. In order to establish the convergence of Theodorsen’s method under reasonably general 
conditions, we employ the integral of Lebesgue. 
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convergence of F,(e*) to F(e*) implies that F,(z) converges to F(z) uniformly for 
| z| <1, and thus the functions f,(z) =ze*"“) converge uniformly for | z| <1 to the 
mapping function f(z). 

In order to prove the convergence of the functions 0,(¢) and 6,’ (¢) we shall derive 
estimates for the differences | On(b) —O(¢) | and | 0, (¢) —0’(¢)| in terms of € and n. 
These differences approach zero as n— ©, and will at the same time permit us to ap- 
praise the degree of accuracy of the mth approximation. 

REMARK. Theodorsen considers the case where the exterior of a circle |¢| =R is 
mapped onto the exterior of a “nearly circular” closed curve T whereby the mapping 
function w= g(f{) is so normalized that lim;.,. w/f =1. This case is immediately reduced 
to the one considered above by means of the transformations w=w~! and z= R/¢. Let 
us suppose that [ is represented by the equation r=r(@) (00 S27), where, for some 
positive b and 0<e<1, b(1+¢)-!S7r(@) $b(1+¢€) and |r’(O)/r(@)| Se. Then the 
function w=f(z) =1/g(¢), where {=R/z, maps the circle | z| <1 onto the interior of 
the nearly circular contour C represented by the equation p=p(6) =1/r(@), where 
@=—0 and p(@) satisfies the conditions (2.1) and (2.2) with a=b~'. For ¢= Re*¥, we 
write arg [g(¢)/¢] =O(W) —y, where arg [g(¢)/¢] is defined as 0 when ¢ = ©. Then, for 
¢=Re’¥ and z=e* where d= —y, 


g(f) 
log - 





log r[O(y)] + iO) — ¥) — log R 


~ 





) 
— log — log R 


VA 


~ 


— log p[6(¢)] — i(0(@) — ) — log R. 
Thus one can form the successive approximations 9,(W) for the function O(y) in the 
same manner as the @,(¢) are formed for 6(¢). Furthermore, 0,(W) = —0,(¢), Y= —@ 
and 0,(~) —O(W) = —(@,(¢) —8(¢)). Hence any bound obtained for | 0.(@) —6(¢)| is 
also a bound for | 0,(y) —O(y)|, and the same remark applies to the derivatives of 
these differences. 

3. Statement of results. We shall prove the following estimates: 

I. If Cis a nearly circular contour, and if 0,(@) and 0(¢) =arg f(e**) are defined by 
(2.6) and (2.4), respectively, then 


2 1/4 
| @n(¢) — 0(¢)| < 2(—*—) emt 212, (3.1) 


1—-—eé 


The bound for | 0n(@) —O(¢) | obtained here approaches zero as n— (since 
0<e<1) and is therefore sufficient to establish the convergence of the functions 
6,() to 0(@). However, a bound which converges to zero more rapidly can be found if 
a further assumption regarding C is made. 

Il. Lf Cis a nearly circular contour and if 


p’(82) p’ (81) 





| | 
| ~ | sel 0.—4|, (3.2) 
} p(2) p(8;) | 
€ being the same as in (2.1), then 
| On(@) — O() | S (24A (m + 1)) Per, (3.3) 


2 
where A =4‘e*. 
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The following result is obtained for the derivatives 0,/ (¢). 
III. If C is a nearly circular contour, if (3.2) holds, and if p(@) =d[p’(@)/p(0) |/dé 


satisfies the condition 


| (02) — p(6:)| < «| 62 — |, (3.4) 
¢ being the same as in (2.1), then 
| 0! (6) — 0 (6) | S V/2mo, (A(m + 1))*ent, (3.5) 
where A =4*e* and 
neteey- wnti¢ela+ wee 3.6) 
For all n, ~ 
on S (1+ €) exp [2e%/xA (1 — €)-*/?], (3.7) 


so that o, 1s bounded if 0<€<1. 
Estimates for the difference | F.(z) — F(z)|, | s| <1, may be obtained from those 


for | 6,(¢) —0(¢)|. For by (2.2), 
| F.(e*) — F(e'*) | S {e2(On1(¢) — 0(¢))? + (n() — 0(6))?} 12, 


and for | z| <1 


| F(z) — F(z) | < max | F.(ei*) — F(e'*) |. 


Thus, for example, in case II we find by use of (3.3) that 

F(z) — F(z) | S 2(Aw(m + 4)". 

Hence, if 0<e<i, the successive approximations F,(z) converge uniformly to 
F(z) =log [f(z)/z] when | z| <1. An analogous statement applies to the derivatives 


d| F,.(re‘*) |/d¢ and | F(re‘*) |/d¢. 
To prove the three theorems I, II, and III, we shall first derive bounds for the 





square means 


2 1 - } 2 72 1 ” , ’ 2 
M, = ~ f (6n() ka 0(¢)) dd, M,, — —f (6,() = ¢ (¢)) d@, 
0 2a 0 


2 


on (3.8) 
y72 ” ” 2 
Mi," = — J (0,'(6) — 0'"()) "de. 


The above results will then be obtained by use of the inequalities (see §4(c)) 


| 0.(¢) — 0()| S (2mM,M,’)'2, —-| 6, (@) — O(¢) | S (29M Ma’)*?. 


The functions f,(z) =ze?*“) map the circle |z| =1 onto closed curves C,. Since 
the functions f,(z) are to be used as approximations to the mapping function f(z), 
it is essential to know that the C, are simple closed curves. This will certainly be the 
case if the C, are star-shaped with respect to the origin. (A closed curve is star-shaped 
with respect to the origin if every ray from the origin intersects the curve in exactly 
one point.) Knowing that C, is star-shaped has the additional advantage’ that 0,(@) 





3 Cf. Theodorsen and Garrick, I.c., pp. 184-185. 
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is then an increasing function of @ and therefore possesses a unique inverse function 
¢=¢,(8). This permits us to form immediately the inverse z=e** of the mapping 
function w=f(z) for w on C,. We examine therefore the question when the C, are 
star-shaped, and obtain the following result: 

IV. If C is a nearly circular contour and if the condition (3.2) is satisfied, then the 
curve C; is star-shaped with respect to the origin if ¢S<(2 log 2)—', C2 if «50.34, C; af 
€<0.31, and C, if «=0.3. For n24 all C,, are star-shaped if €30.295. 

This result is derived by examining the values of ¢ for which | 0,/ (¢) — 1| <1, so 
that 6,’ (¢) 20 and 6,(@) is therefore monotone increasing. For large values of (n= 4) 
a more favorable estimate for € may be obtained by making use of (3.5) and of a lower 
bound for @’(¢) which is given in §9(d). 

4. Proof of I. (a) EstimaTE oF M,. Let F(e**) and F,(e*) be the functions in 
(2.4) and (2.7). Because of the representations of 0(¢) —¢ and 6,(¢) —¢ by means of 
the integrals (2.5) and (2.6), respectively, we have 


f "(6() — db = 0, f “(Ga(d) — 6)d6 = 0. (4.1) 


We now apply the following well known theorem If the function g(@) is real- 
valued, periodic (period 2m), and (g(@))* is integrable (in the sense of Lebesgue) 0S S2r, 
and if 2(p) is a conjugate function of g(p) (then surely existing), then 


1 Qn 1 Qe 
— J le) Fe + at = — foto }rae + 6° (4.2) 


0 


1 2r 1 2x 
a=—f[ ‘sede, b=—f goas. 
2a 0 2r 0 


Applying this with g(¢)+i2() = F,(e*) — F(e*) and observing that B=0 (be- 
cause of (4.1)), we obtain 


2 1 2r 2 1 Qe 2 
M, = — f (0.(¢) — 0($)) de S — f {log p[4.-1(¢)] — log p[o(d) }}ag. (4.3) 
2r 0 2n 0 


where 


By hypothesis (2.2), 
| log p[@n-1(¢)] — log p[0(¢)]| < €| A.-1(6) — 0) |, 


and therefore 


2 2 1 ” 2 2_.2 
M,, ~ re (6,-1() -) 6()) do 4 € Ma~ts 
2 0 


or 
M,, =< eM,-1, M, s M ce". 


For »=0 we obtain from (4.2) by use of (2.1), 


2 1 2e ° 1 2r 6 2 
M, = —{ (Q(¢) — ¢) do S —{ (10g ) dé Ss ¢. 
2 0 2r 0 a 


4 See, for example, A. Zygmund, Trigonometric series, Warsaw, 1935, p. 76 (Eq. (4)). 
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Thus we have proved that if p(0) satisfies hypotheses (2.1) and (2.2), 
M, & e"™*', (4.4) 


(b). Estimate or M,). It follows from §9(b) and §9(c) that F,(e**) and F(e**) are 
absolutely continuous and that {d[F,(e*) ]/d@] }2 and {d[F(e*) ]/do }2 are integrable. 
Furthermore, because of the absolute continuity of F,(e‘*) — F(e*), the imaginary 
part of the derivative d{ F,(e**) — F(e*) } /d is a conjugate function of the real part. 
Finally, 


2r d outer Qn d 
J Grtende = (ren li = 0, f  <Fatenae = 0. 
0 0 do 
Hence, applying (4.2) with g(¢) +i2(¢) =d[F,(e*) — F(e*) |/db, we obtain’ 


2 1 - , , 2 
3 f (01) — 0'(6))"de 
2rd o 





SF at A , p’ ™ 
=f 4% -wor|otsc) - = ore e'ce} as. (4.5) 
2rd o p p 
By (2.2) we have (omitting the argument ¢ in the integrands) 
sx f" (0... + 0” )do. (4.6) 
As is shown in §9(c), 
1 2 
—f 0d S (4.7) 
rd 9 1—é 


Furthermore, applying (4.2) with g(¢)+ig(¢) =d[F,(e**) ]/dd, we obtain by (2.2), 


1 
fe — 1) “do = —f (e (6,110, ) dp Ss Se f° 6, dé, 
2r 

1 2r 

—{f aap —2f ‘ade + 2b sof 0. de. 

2a 0 0 


Since Sr"0. dd = 27, we find that 


1 Qe y 1 Qe . 
me f i or f 0. de. 
2x 0 2x 0 


or 


If we set 


we have m,?S51+m,2_,e, and therefore 
2 2 4 Qn 2 
mM, Site te +--- +e m. 


Since 


5 The notation “[0] or (p’/p) [a] means p’(6)/p(@). 
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1 


‘ 2 
Mo = - : do = 1, 
we obtain 
mS (1—-e¢) . (4.8) | 
Thus by (4.6), (4.7), and (4.8), 
we be HT ~ 9). (4.9) 


(c) EstiMaTE oF |6,(¢) —@(¢)|. To complete the proof we now apply the follow- 
ing theorem: Jf g(¢) is a real-valued, absolutely continuous and periodic function 
(period 21) and tf (g'(p))? ts integrable, then for any do, 


[2(¢) ]? — [g(¢o) }? S 2nMM’, (4.10) 


where 


as f [e(d) |", 9 M2 = — f [e’(¢) }?a¢. 
2rd 2rJ 
The factor 27 is the “best possible” constant; it cannot be replaced by a smaller one.® 

Let g(¢) =9,.(¢) —9(¢). Since then S"g(0)do =0, there exists a value ¢o such that 

g(¢o) =0. Hence 
| @n(@) — 0(¢) | S (24M ,My)"/?. 
Using (4.4) and (4.9), we find (3.1). 

5. Proof of II. (a) EstimaTE oF M,;. Under the present hypotheses an estimate 
for M,; sharper than (4.9) may be obtained. We shall prove that if p(@) satisfies (2.1), 
(2.2) and (3.2), then 
M,' <A(n+ te, = (A = 4¢e*). (5.1) 


Using the relation (4.5) we obtain for n21, 


= fi T(E r eo yy eo oe ny. ‘ 
M, it ; [6,1] : [eo] jo’ + : [0,1 ]}(O/_1 — 6’) | do 


sf" (Eton -£ yoru} 


ff , < oe Pe 2 \ 1/2 
+ lon . (6, —1 6’) , [0,1] do ’ 


6 To prove (4.10), we note first that for OS ¢S27, OS goS2z, 


¢ o-2" 
g2(¢) — g2(go) = 2f g(t)g'(t)dt = 2f g(t)g’(t)dt. ) 
° ?0 do 
Since 
¢ | po-or ¢ or 
gg’| dt) + f gg’| dt -f | 68’| a= f | ge’| dt, 
“"¢ 1 I$, o—2e 0 
one of the two integrals in (*) does not exceed 4/°" | pp’| dt. Hence, by the inequality of Schwarz, 


or 
[e(¢) ]? — [g(¢o) }? S f | ge’| dt < 2nMM’. 


Applying (4.10) with g(¢) =cos" ¢(¢0=42) and letting n— ~, we see that the constant 27 cannot be re- 
Ppl) g 


placed by a smaller one. 
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by Minkowski’s inequality. Under the present hypotheses we have by §9(d), 
0<6(¢) SA = 4e*, 
Hence, by (3.2) and (2.2), 


1 Qe 1/2 
M, s €A -S (0,1 = aad + eM i_, = €(A M,-1 + M,_,), 
TH 0 


and, therefore, by (4.4), 
M,,’ =< e(Ae” + M,.1), (n = 1). (5.2) 


For n=0, we have, using (4.2) and (2.2), 


1 Qe 1 Qn p’ 2 A Qe 
My? = f (0 — 1)%dp = —{ — [ele )des eX f dp = A. 
2rJo 2rd o p 2rd o 


This inequality proves (5.1) for »=0. For n21, (5.1) is easily seen to be true by 
induction. Assuming that it holds for some n20, we obtain by use of (5.2), 


Masi S €(Ae"t! + A(n + 1)e"*’) = A(n + 2)e"*?, 
i.e., (5.1) is also true for n+1. 

(b). Estimate oF |6,(¢)—6(¢)|. Applying (4.10), (4.4) and (5.1), we find 

| On(@) — 0(¢) | S (2mM,M,/)"!? S (2eA(m + 1))*/2e"+1, 

6. Proof of III. (a). SOME PROPERTIES OF THE FUNCTIONS F(e**) AND F,(e**). Be- 
cause of the hypothesis (3.4), F(e**) has a continuous second derivative for? 0 <@ S2r. 
The same is true for all F,(e**) as is shown in §9(e). Differentiating F(e**) and F,(e**) 
twice with respect to ¢, we obtain 


U 





dF’ °F p 

— = — [0|0’ + ie’ — 1), — = p|0}0? + — |6)0” + 16”, 
Pate [2] ( ) ie p|é] : [6] 

dF, p’ ’ aii n 72 p’ ” ont? 
db - y [0n—1 ]On—1 + i(6, = 1), dep? = P[On—1]0n—1 +> a [4n—1 JOn—1 a 18n. 


The present proof is similar to that of (II) and we estimate first 


1 2r 1/2 
M.’ = \ f (6! — ayragh 
T/ 0 


Mi! < A%(n + 1)%one"*1, (6.1) 


We prove that, forn2=1, 


where A =4*e* and a, is defined in (3.6). 
(b). PROOF OF THE INEQUALITY (6.1). Since 


Qr a? 
f <" (ei) — Fle'*))de = 0, 
o dg? 


Trans. Amer. Math. Soc. 38, 326 (Theorem III), (1935). 
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we have, applying (4.2), 
72 


tT 2 1 ” 72 72 
(Me) = = f {ole = aloe” + plo." — 6 } 
+ (* “eo (o])0” +— Ble - a” do. (6.2) 
p p p 


Because of (3.2), 
| p(0)| Se. (6.3) 


Using (9.3), (3.4), (6.3), (3.2), and (2.2), we find that 


us 


. 2 2r . . 
(Maw) S < f {A°|0, -—0|+|o, —0” | +|0”||o,—0| +| 0’ — 0’ |} de. 
0 


If M, and M,’ are defined as in (3.8), we have by Minkowski’s inequality: 


” 2 1 ae 2 se 
Ae {4 M, + ra (6, — 6 ae} 
Qn 0 


1 - rn2 2 “_ ” , 


lr 


1 Qe 1/2 a 
M3' = if (o”yraa S A*/*%/2e, (6.5) 
TH 0 


| @n(@) — O() | S (2A(m + 1)) Per}, 


Since by (9.4), 


and by (3.3), 


we have 


‘1 Qe 1/2 a 
a (0, — aarrags < (2"A(n + 1))/2A3/24/2 ent? = 24%™+%/a(n + 1). (6.6) 
0 


2r 
Next, applying the theorem of §4(c) with g(¢) =8, (6) —80’(@), we obtain 
(0/f — 6)? 5 2xM, M,Z’, 
and taking the square root and using (9.3), we have 
10. +6 | S24 +/2eMi My’. 


Hence 


1 2r 1/2 wees fe oe ee. Se 
{ f (6/2 — vsyagh < Mi! (2A + 4/2eM/M") = 2AM! + Mis/2eMi MS. 
0 


2r 
Applying the inequality® M,) S\/M,M,!' to the factor M, of the square root we 
find that 
8 If we set g(¢) =0,(¢) —0(¢), we have by integration by parts 
en - ¢* | 
f wonde=| eel - f lore"@ree | 
0 ' 0 
wal 2 a 1/2 
sf" soe") dos ff lobes f° er onrel 
0 0 0 


This proves the inequality of the text. 
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1 2r 1/2 
{— f (0/2 — yah < 2AM/ + Mi's/2nMM. (6.7) 
To 


Thus we obtain from (6.4) using (4.4), (6.7), (6.6) and (5.1), 


Mis S €{ A%e"+) + 242m + 1)e"*! + 2A %e"+24/e(m + 1) 
+ (1+ e+4\/2eA(n + 1))M,}, 


and therefore 


7 


ee ities: - a 
Mis ed + 2(n + 1) + 2er\/a(n + 1) + (1 + *+'4/2eA(n $1) \. (6.8) 


Arenti 





Assuming now that (6.1) is true for some 22, we see from this inequality that 
(6.1) also holds for +1. For, if we substitute in (6.8) for M,!’ the right-hand side of 
(6.1), we find that 


Miia S Aten*?{ 1 + 2(m + 1) + 2ev/a(n F 1) + (1 + ett 2A (n $ 1))(n + 1)%0n}. 
For n=2, 
1+ 2(m+ 1) + 2e/e(n +1) < [14+ 2%n4+ )](1 +6 < [1+ 2m 4+ 1) Jon, 
and therefore 
Mis S A%e"*Gngi(1 + 2(m + 1) + (nm + 1)2) = A2(m + 2)?onys€"*?. 
To complete the induction we show that (6.1) holds for »=1 and n=2. From 


(6.2) with »=0, we find that 


1 Qe : . 
M,' = ~ f | (oo) — p[0(¢)])6’? + p(o)(1 — 0%) — * (oo)]e" io} 


1 2r 1/2 
ae My + (~ f (i — 46) + us 
Qn 0 


by Minkowski’s inequality and (3.4), (6.3), and (2.2). Applying (4.4), (9.3), (5.1), and 
observing that by (9.4) Mg’ SA*e(1+€), we find that 


1/2 


IIA 


1+A 
Mi’ < &(A2?+ (1+ A)A + A140) = Aré(2 $ ate + ‘). 


Since A 21, (1+A)/A S2, and therefore 
Mi’ S A%e?(4+ €) < 4A7e?(1 + €). (6.9) 


To prove (6.1) for »=2, we apply (6.8) with m=1 and My’ replaced by A*e?(4+€) 
(see (6.9)), to obtain 


Mi! S A%e[5 + 2ey/2e + (1 + 2€%/wA)(4 + ©]. 
Since 2\/27 <6 and 1+2¢€\/rA >1, 


5 + 2ey/2e + (1 + 2e%/eA)(4 + ©) < (SH 66+ 4+ €)(1 + 2€%/xA) 
< 91 + 6)(1 + 2€V/xA) = 3 on, 
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and, therefore 
M3’ S 3%02A7e?. 
(c). EstTiMATE oF | 6,’ (¢)—0’(@)|. Applying the theorem of §4(c) with g(@) 
= 6, (6) —0’(¢), we obtain from (5.1) and (6.1) 
| 0x, () — 0'(¢) | S V/ 2m, (A(m + 1))?/ent?, 
(d). Proor or (3.7). To estimate o, we first note that 


I] (1 + *\/27Ak) S exp | ver > vi, 


k=2 


Now 


n n n n 1/2 
te e*n/k = € > e(F-1)/2(4 ReF-DI2) < f = eF-1 >» beh : 
k=2 k=2 k=2 


by the inequality of Schwarz. Hence 


n : 1 1/2 2, /9 
Devise} : ( -1)} Poe... 
par 1—e \(1 —.)? (1 — «)3/? 


We find therefore that ¢,<(1+¢) exp [2/xA e?(1—6)-*/?]. 

7. Anintegral representation for 0,/ (6). We shall discuss now the conditions under 
which the images C,, of the unit circle by means of the functions w=f,(z) =ze**™ are 
star-shaped. For this purpose we shall first establish the following representation for 
6, (6). If C is a nearly circular contour and if the function p(@) which represents C satis- 
fies hypothesis (3.2), then the derivative 6, (@) of 0.(@) ts continuous and 











1 ore r 4 t— 
6{(@) -1=-—- —{ {f (i) — a co} cot ° it (7.1) 
2rJ gx Vp p 2 
1 re (,’ p’ i—@ 
6/(¢) -1=—-— + (é,-1(é)] — — .-s@|} 6,/_,(t) cot dt 
2rJy \p p 
p" p" , ” 
ed ae [0,-1(¢) | _— [0,-2(¢) ]6,._2(¢) (n = 2). (7 : 2) 
p p 


Proor. The integrand of (7.1) is continuous in both the variables ¢, ¢ except pos- 
sibly for =, and is bounded because of (3.2). Hence the integral (7.1) is a continuous 
function of ¢. Since this integral represents the conjugate function of (p’/p)[@] for 
which the integral over the interval (0, 27) is zero, it is equal to 0 (¢) —1, at least for 
almost all ¢, and, because of the continuity, for all ¢. This proves (7.1). 

Let us suppose it were proved that 6/(@) is a continuous function when 
k=1, 2,-+-,m (m21). We then show that the formula (7.2) holds with » replaced 
by +1, and that 6,/,:(¢) is continuous. This will then prove the representation (7.2) 
and the continuity of 6,’ (@) for all n. 

Since Fy4:(e**) =log p[6n(b) | +2(On41(¢) —) is absolutely continuous (see §9(b)) 
it follows that 0,’4:(@) —1 is conjugate to (p’/p) [6,(¢) ]@,/ (), and we have, for almost 


all ¢, 
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T=o+t 


1 * (p’ t 
6x41(@) a 1 alien —f {= [6,(7) Jorn} cot — dt, 
2rJ o p 2 


t=o—t 


the integral being convergent in the sense that lims.o {jf exists. We write 


1 ¥ (0! ; t 
6/i:(@) -1 = — —{ {= [0,.(@ + #)] — oa (oxo) 6, (@ + t) cot — dt 
2rd o p p 2 


1 ® U , t 
a ‘+ [0.(@ — )] — dl once) 6,, (@ — t) cot — dt 
24 p p 2 


p’ 1 aa P , » & 
- tools f [0 (@ +) — 0: (6 — )} cot— at 


Because of (3.2) and the continuity of 6,’ (@), the first two integrals represent continu- 
ous functions of ¢. The third integral (without the factor —(p’/p) [@,(@) ]) is equal to 
(p’/p) |0,’ -1() |0,f _1(@), since 6,/ (@) —1 is conjugate to this function. Introducing the 
variable r=¢++¢ in the first integral and tr =@—1 in the second, we obtain 


’ = ce LJ iain e m. ‘A ’ i—@¢ 
6/.4(¢) -—1= = + [6,(7) | - onc) 6 (r) cot 





dr 


—* [6,(6)] ~ [0,-1(6) ]6x’-1(4). 
p p 


The right-hand side of this equation represents a continuous function of ¢. Hence 
6,/,1(@) may be defined as a continuous function for all ¢, and therefore 0,4:() has a 
continuous derivative for all ¢. This completes the proof. 

8. Conditions under which C,, is star-shaped. Proof of IV. The curve C, is star- 
shaped if 0,7, (@)20. By (7.1) and (3.2) 








ire|/ ae |¢—@ 
a) -1ls—f — (t) — — (¢)| cot | ——| dt 
2rJ 4. | p p | { 2 
7 ott i- ’ 
< “s (t — @) cot dt = 2¢ log 2. 
Tv o—4r 


Thus, if 2¢€ log 2 <1 or eS(2 log 2)—, then 6/ (¢) 20 and C, is star-shaped. 
Let us suppose it were proved that 6,’ (¢) 20 for some 21, provided € does not 
exceed some value €)<1. Then we examine 6,4:(¢). By (7.2), (3.2), and (2.2), 


i- 
oute 








€ ot+r 
|g? = Bal = ’ 
|ofss(6) - 1] $= J Gd ~ 4(6))8 ( cot /1(6)|. (8.1) 


TT o—r 


It is to be noted that 0,(t)—8@,(¢) has the same sign as t—@ since 6, (t) 20. We find 
by integration by parts that 


» 1 fee 7 t—¢ 1 fr" On(t) — On) ) 
n ae On = On 6, at => ad d 
a J og Cee? — GhEI O Cot —, oi... Vana 


ns = Me —t— [6.(¢) -—¢]) +¢- Va 
2a J os 2 sin $(t — ) 
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Hence by Minkowski’s inequality, 


i OD ll te —— [2.(6) — a 
2rd 4s 2 sin 3(t — ¢) ) 


1 ott i—¢ ae 
+e. Gawcal 
2xaJ 4. \2 sin $(t — >) 


Integrating by parts, we find that 


1 o+r i - ¢ 2 1 o+r t{— ry 
— ; —*—) dt = —{ (¢ — ¢) cot - dt = 2 log 2 = &’. 
2e/J4-- \2 sin $(t — ¢) 2a J gx 2 











Furthermore, by the theorem of §9(f), 


1 i —+t— [6,.(¢) - oN a ae nie p[4n-1(t)] — log cee NNN a 
2 sin 4(t — ¢) ~~ 2 sin 3(¢ — ¢) , 








2a J gx 2 


By (2 2), the right-hand side of this equation is 
2 1 “CS — On-1(¢) ): 2 2 
— dt 


= € Mn-1. 
2 sin 3(¢t — ¢) 





IIA 


€ 
2r o-—F 


Hence 
My S €My-1 + C. 
Since mo =c, we have 
1 an e™tl 


mscitetet+-:---+e)=c 


1—e 


Hence, by (8.1), 


lA 


1 a ettl 2 2 
2 j=] log 2+ | Oi1(¢)|. (8.2) 
€ 


Of41(¢) -—1| S em, + e | 6,1-1(¢) | “ie 





Applying (8.2) with n=1, we find since 6¢ (¢) =1 that, 
| 02 (@) — 1| S 2e(1 + €)* log 2+ e, (8.3) 


and this will be less than 1 if €$0.34. 
For n=2 we find, since 6/ (6) $1+2e log 2 and @/ (@) >0 for €<(2 log 2), that 


| of (@) — 1| S 2e(1 + € + €*)* log 2 + e(1 + 2c log 2). 


This expression will be less than 1, if e€$0.31 (<(2 log 2)~). 
By (8.3), | 6, (¢)| $1.7927 if €=0.30. Hence, applying (8.2) for »=3 and using 
this estimate for 67 (¢), we find that | ai (o) — 1| 31, if eS0.3. 
Assuming that, for some »21, 0<8,/_1(¢) $2, we see from (8.2) that 
2 log 2 
ame 


if ¢<0.295. Since for «<0.295 this assumption is certainly satisfied for n=1 and n=2, 
it follows that for all m |6,/4:(¢)—1| <1 if eS0.295. 


e+ 2e2< 1 





| Ox1(@) —1| S 
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REMARK. For large values of m the bound for ¢ can be improved by use of Theorem 
IV and the left hand ineqyality in (9.3). 

By Theorem IV, |6,/(¢)—6’(¢)| S$/2x0,(A(m+1))*/2e"!, and by (9.3), 6’) 
=A-(1+¢)-/*. For any given fixed n, €) can be chosen so that \/20n(A (n+1))*/2eG*" 
<A-(1+¢)—/*. Then, for all € Seo, 0,/ (¢) 20’(¢) —A“(1+€2)—220. 

9. Auxiliary theorems. This section contains the proofs of some of the auxiliary 
results cited in the text. 

(a). UNIQUENESS OF THE SOLUTION OF THEODORSEN’S INTEGRAL EQuatTion. Jf C 
is a nearly circular contour, the integral equation (2.5) has at most one continuous solu- 


tion. 


Let us suppose that it had two such solutions, 0:(@) and 62(@). Since 
Qe 2r 
fe -ode=0, fee) - #49 =0, 
0 0 


it follows by use of the theorem cited in §4(a) that 
1 2r 1 2r 
M? = — f (0:(@) — 62(¢))*dp < — f {log p[4:()] — log p[62(4)) |} *de. 
24/0 24 Jo 


By (2.2), 
| log p[0:(¢)] — log p[62(p)]| Se] A(d) — 62()|, 


so that we have M*Se?M*. Since 0<¢€<1, M=0 and hence 6:(¢) =62(¢). 


(b). A PROPERTY OF THE FUNCTIONS 6,(@). If C is a nearly circular contour, then the 
functions 0,() defined by (2.6) are absolutely continuous, and (6,' (p))* are integrable 
(in the sense of Lebesgue) for 0S¢S2r. 

This is clearly true when n=0. We suppose that this statement were proved for 
some 20. Since log p(@) has bounded difference quotients (by (2.2)) and @,(@) is ab- 
solutely continuous, it follows that log p[6,(@) | also is absolutely continuous. Further- 
more, because of the inequality 


(2 lo(6) as )) 5 C@L (6) 


it follows that the integral 


f as (= [6.(¢) ]6." (@)) ae 


exists. Hence, the conjugate function of log p[@,(¢)], namely 6,4:1(¢) —@, exists and 
is absolutely continuous and the integral 3" (0.41(¢) —1)*d@ exists.® 


(c). A PROPERTY OF 0(¢). If C is a nearly circular contour, then 0() =arg f(e*) (de- 


® We are using here the following theorem: if g(¢) is an absolutely continuous and periodic function 
(period 27) for OS ¢S2- and if [g’(¢) FP is integrable, then the conjugate function g(¢) is absolutely con- 
tinuous and (’(¢))* is integrable. (See, for example, W. Seidel, Uber die Rinderzuordnung bei konformen 


Abbildungen, Math. Annalen 104, 223 (1931). 








26 S. E. WARSCHAWSKI [Vol. III, No. 1 


fined by (2.4)) is absolutely continuous and (0'())? is integrable (in the sense of Lebesgue) 


and 


2r 


1 
= 0°(¢)do S eee 


2nrvJ o —é 


(9.1) 





ProorF. Since the curve C is rectifiable, the function F(e‘*) is absolutely continu- 
ous.’ Hence 


d . p’ 
—F(e*) — i = — [0(¢)] 6’'(¢) + 16’(¢) 
do p 


exists almost everywhere for 0 $¢ S27, and is integrable. Furthermore, the function 
d[ F(z) |/d¢ —i=u(z)+i0(2z), z=re*, may be represented by the Poisson Integral in 
the unit circle, 


2 1-—r 


1 2r 
u(z) + iv(z) = —f } u(et*) + iv(e**)} er dt. (9.2) 
0 


2r r? — 2r cos (t — ¢) 





For almost all ¢(0S@ S27), 


r—1 


lim u(re’*) = .. [0(¢) |0’(¢) = u(e**), lim v(re'*) = 6’(o) = v(e'*). 
p 


r—1 
Since C is star-shaped, 0’(@) 20, and we have by (2.2), 
v(e*) + u(et*) = 0’()(1 — «) 2 Oz 


Because of the representation (9.2) we conclude that v(z)-+u(z) 20 and v(z) —u(z) 20 
for | z| <1. Hehce v?(z) —u?(z) 20 for | z| <1. Now 


II 
— 
: 


1 Qn 
—{ (v?(re‘*) — u*(re*))dp 
2rd o 


Hence, taking the limit as r—1, we obtain by Fatou’s lemma 


1 2r p’ 2 
— 6’ 1—{— d 
al @)| ( (oe!) | 6 


2 1 
f (6'($))*dd S 


11 
2rd o 1 —e? 


IIA 


1, 


and by (2.2), 





This proves that (6’(@))? is integrable and that (9.1) holds. 


(d). AN ESTIMATE FOR 0’(¢) AND 0’’(@). If C ts a nearly circular contour and tf in 
addition (3.2) is satisfied, then 








J S64) $A = Ate", (9.3) 
A J/1 a €? 


1 Qe 1/2 - 
\— f | 0"(¢) } ‘ah < A*/%e min (1 + €; 4/2). (9.4) 
2r J 9 

1 This follows from a theorem of F.and M. Riesz, Uber die Randwerte einer analytischen Funktion, 
Comptes Rendus du Quatriéme Congrés des Mathématiciens Scandinaves 4 Stockholm (1916) pp. 27-44. 
See also, F. Riesz, Math. Zeitschrift, 18, 95 (1923). 


1945] CONFORMAL MAPPING OF NEARLY CIRCULAR REGIONS 27 


The proofs of these inequalities are contained in a paper to be published elsewhere. 


(e). A PROPERTY OF THE FUNCTIONS F,(@). If C is a nearly circular contour for 
which (3.2) and (3.4) are satisfied, then the functions F,,(e**) have continuous second de- 
rivatives which satisfy a Holder condition with any fixed exponent a, 0<a<1. 


The proof may easily be given by induction. Since log p(@) and 6;(@) —@ are con- 
jugate functions and since the second derivative of log p(@) satisfies the Lipschitz 
condition (3.4), it follows from a theorem of I. Privaloff," that 6/ (@) and 6{’(@) exist 
and that 6/’ (@) satisfies a Hélder condition with any fixed exponent a, 0<a<1. Let 
us suppose now, that it had been shown that 0,'’ (¢) exists and satisfies a Hélder con- 
dition with any fixed exponent a, 0<a<1. Then log p[8.(¢) | has continuous first and 
second derivatives, (p’/p)[6.(@)]0./ (6) and b[8n(¢) 10,°() +(p’/p) [0.(@) |0,’’ (@), re- 
spectively, and the latter satisfies a Hélder condition with any exponent a, 0<a<1, 
(because of (3.4) and (3.2)). Hence, again by Privaloff’s theorem, the conjugate func- 
tion 0,4:1(¢) —@ possesses a second derivative 6,',,(¢) which satisfies such a Hélder 
condition. This completes the proof. 

(f). A THEOREM ON CONJUGATE FUNCTIONS. Let us suppose that u(t) is a periodic 


function (period 2m) possessing a continuous derivative for OStS2n. Let v(t) be con- 
jugate to u(t) and let us suppose that v(t) also possesses a continuous derivative for 


0<t<2r. Then for every 6, 


Qr i) -— 6 2 2r t)— a] 2 
f (= Vat = f (=< no) at 
0 sin 3(t — @) n sin 3(t — 6) 
Proor. Let G(z) = U(z)+iV(z) denote the function which is analytic for | s| <1 


and assumes the boundary values g(t) =u(t) +40(t) — (u(@) +-10(8)) for z=e*. Then the 
real part of [G(z)]* may be represented by the Poisson integral (2 =re‘) 


1-—r 


1 Qe 
este fered an”. F a = 2 
[U(z) ]? — [V(2)] f { (u(t) — u(@))* — ((2) — 0(6))?} ish-beneoa dt. 


Oe 





For z=re*®, 








[U@Pr- Vert f" (u(t) — u(0))* = (0(H) = 00)? 


ee "fe (1 — r)? + 49 sin? 3(¢ — 8) 


As is easily seen, the limit of this integral as r—1 is 


1 ¢-** (u(t) — u(6))? — (v(t) — v(8))? a 


— - (9.5) 
2rd o 4 sin? 3(¢ — @) 





By the mean value theorem (since U(e**) = V(e*) =0) 
1 |. Privaloff, Sur les fonctions conjuguées, Bull. Soc. Math. France 44, 100-103 (1916); or Zygmund, 


l.c. p. 156. Privaloff’s Theorem states: if g(@) is periodic (period 27) and satisfies a Hélder condition with 
the exponent a, 0<a<1, for all ¢, then any conjugate function of g(@) satisfies such a condition. 
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Ui)? —-[V@pPpa,. . go * at 
or — = — [U%(pe) — V%(pe*) | nz 
1-—r,r Op 
r) ) ; 
a {Uloei — U(pe*) — V(pe*”) — Vise) »(r<#% <1). (9.5) 
Op Op 


p=r 

Since g’(#) exists and is continuous, 

oU _ OV aoe : : 

— t+ §— = e*G"(pe*) — — ig'(0) 

Op Op 
as p—1. Thus lim,.: 0[U(pe*)]/dp and lim, 0[V(pe*)]/dp exist. Furthermore, 
lim,.1 U(pe**) =lim,.1 V(pe**) =0. Hence, the limit as r—1 of (9.6) is zero and there- 
fore the integral (9.5) is zero. This proves the theorem. 
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ON ROTATIONAL GAS FLOWS* 


BY 


ANDREW VAZSONYI 
Elliott Co., Jeanette, Pa. 


Introduction. The main body of the science of aerodynamics is based on the classi- 
cal theory of frictionless, incompressible, irrotational fluids. Recently airplanes have 
attained such high velocities that this fluid model has proved to be too restricted and 
interest has centered on the irrotational motion of frictionless, compressible fluids. 
By the term “compressible fluid” one generally means a fluid for which the density 
p and pressure p are connected by the isentropic relation pp-’ =const. However, the 
student of aerodynamics is frequently interested in supersonic phenomena and be- 
cause of the possible occurrence of shock waves, such flows cannot be described, in 
general, by isentropic, irrotational flows. Accordingly, it becomes necessary to study 
the motion of gases under less restricted conditions. 

Let us call a fluid barotropic when there is a unique functional relationship be- 
tween the pressure and the density of the fluid. The most important examples are the 
incompressible fluid where the same constant density belongs to each pressure and 
the isentropic fluid where the relation pp~7¥=const. holds. The dynamics of friction- 
less barotropic fluids is based on a theorem due to Lagrange. If a fluid particle is 
irrotational at one moment, it will remain so for all subsequent time. One can generally 
assume in aerodynamics that the air starts from rest. The dynamics of the flow can 
then be summed up in the single statement that the motion is irrotational. It follows 
that the velocity distribution admits a potential, and the comparative mathematical 
simplicity of the dynamics of frictionless barotropic fluids follows from this fact. 

Classical fluid dynamics deals almost exclusively with the theory of frictionless 
barotropic fluids. To find an example of frictionless non-barotropic fluids, we turn to 
the theory of the propagation of waves. When Newton developed his theory of sound- 
waves, he assumed that the motion of air was isothermal. Later his theory was super- 
seded by a better one which assumes isentropic motion. Thus, both theories assumed 
‘that the transmitting medium was barotropic. The mathematical theory of one- 
dimensional large disturbances, a much more difficult problem, was developed first by 
Riemann, who again assumed that the flow is isentropic. However, the isentropic 
theory of shock-waves turns out to be fallacious because it can be shown to violate 
the law of conservation of energy. When shock-waves are considered the fluid model must 
be extended to include non-barotropic fluids. 

In this connection, let us draw attention to the thermodynamical aspect of the 
general theory of compressible fluids. In the case of a three-dimensional flow there are 
six unknowns: three velocity components, pressure, density and temperature. The 
laws of conservation of matter and momentum together with the equation of state 
yield only five equations. To get the missing sixth equation the law of conservation 
of energy, i.e., the first law of thermodynamics, must be used. Flows will be isentropic 
only when as a consequence of these laws the entropy turns out to be a constant. 


* Received Oct. 2, 1944. 
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Although the theory of one-dimensional shock-waves requires a non-barotropic 
fluid model, this fluid model is a very special one. Even if there is an increase of 
entropy across shock-waves, the flow remains isentropic between shock-waves. More- 
over, a one-dimensional fluid motion is always irrotational But when one turns to 
two or three-dimensional shock-waves, the situation becomes quite different. In this 
case Hadamard! was the first to point out in 1903, that vortices are generated suddenly 
by shock-waves and, in general, the flow becomes non-barotropic after shock-waves. 

Hadamard determined the sudden change of circulation across a shock-wave but 
was not interested in the circulation variations occurring in the fluid behind shock- 
waves. A general circulation theorem for frictionless barotropic fluids was established 
by Bjerknes? in 1900, for the purposes of his dynamical theory of meteorology. The 
motion of air masses originating from non-homogeneous conditions is clearly a phe- 
nomenon requiring a non-barotropic fluid model. 

Crocco,’ in 1937, again took up the question of the motion of frictionless fluids 
behind shock-waves. By restricting himself to the steady state he discovered a very 
useful theorem. Recently, this theorem was generalized by the author of the present 
paper.‘ 

So far, we have spoken only about frictionless fluids. There are problems with 
respect to the flow of gases where viscosity cannot be neglected. We mention, for 
instance, the boundary layer theory and the behavior of a gas within shock-waves. 
It appears probable that when considering the viscous flow of gases, the conductivity 
of the gas cannot be neglected in general. Variations in viscosity might have impor- 
tance also. No general theorems are available for such flows and we shall have to be 
content with presenting the fundamental differential equations governing these phe- 
nomena. Any investigation with respect to the flow of gases must be based on these 
equations. While in the case of frictionless flows some general consequences of the 
fundamental equations are available, in the case of viscous flows we must start the 
investigation of each problem by examining the fundamental equations anew. 

Lagrange’s theorem plays a fundamental role in our concepts about fluid dynam- 
ics. Its validity is restricted, however. The art of aeronautics is now at a point where 
we have to extend our fluid model and thus modify some of our basic concepts. We 
must accept for instance the fact that vortices can be generated in the midst of a 
frictionless fluid. Whether this extended fluid model will be able to account for all 
the phenomena which we may wish to consider, only the future can tell. 


I. THE FUNDAMENTAL EQUATIONS 


1, Continuity equation. From the law of conservation of matter it can be proved 
that 
Op 
div pq = ——> 
Ot 
1 J. Hadamard, Sur les tourbillons produit par les ondes de choc, Note III, in Legons sur la propagation 
des ondes, A. Hermann, Paris, 1903, p. 362. 
2 V. Bjerknes, Das dynamische Princip der Circulationsbewegungen in der Atmosphare, Meteorologische 


Zeitschrift, 17, 97-106 (1900). 
3 L. Crocco, Eine neue Stromfunktion fiir die Erforschung der Bewegung der Gase mit Rotation, Zeit- 


schrift f. Ang. Math. und Mech., 17, 1-7 (1937). 
* A. Vazsonyi, On two-dimensional rotational gas flows, Bull. of the American Mathematical Society, 


50, 188 (1944), 
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where q is the velocity vector. Another useful form of the continuity equation is 
given by 
divq=——- (C) 


2. The Navier-Stokes equation. From the law of conservation of momentum it 
can be proved that the equation of motion is given by 


dq 1 A M . 
— = — — grad p + — Aq + — grad div q, (M) 
dt p p 3p 


where it is assumed that the viscosity yu is constant. 
It will be useful to derive certain other forms of this equation. The specific en- 


thalpy h of a gas is defined by 
h=U-+ po", (2.1) 
where U denotes the specific internal energy. The specific entropy s is defined by 
Tds = dU + pd(p~) (2.2) 
where T denotes the absolute temperature. From Eqs. (2.1) and (2.2) it follows that 
Tds = dh — p-dp. (2.3) 
Using vector notation and considering only spacial variations, we may then write 
T grad s = grad h — p~ grad p. (2.3’) 
From this last equation and the equation of motion we find that 


dq Kh fod . , 
— = T grad s — grad h + — Aq + — grad div q. (M’) 
dt p 3p 

Another useful form of the equation of motion can be obtained by using the 
stagnation enthalpy 
ho = h + 3¢? (2.4) 


and the identity 


dq oq 
— = — + grad (}q7) —q X@ (@ = curl q) (2.5) 
dt at 
together with the equation of motion (M’). Thus one obtains 
oq i K . ” 
5 7 1X @ = — Brad ho + T grad s + — Aq + — grad div q. (M”) 
6 p p 


A fourth useful form of the equation of motion can be obtained by introducing the 
rate of change of the stagnation enthalpy. Differentiating Eq. (2.4) with respect to ?, 
we obtain 


dho dh 1 dq? 
meee a Sed mae ee (2.6) 


eo #: 2.@ 


From Eq. (2.3) it then follows that 
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dho ds 1 dp dq _ds 1 Op 1 dq Ms 
Swpelyt Vey Sys Sel gg eyes en 
dt dt p dt dt dt p at p dt 


~ 


Using the continuity equation (C’), we may write this as follows: 


dho ds 1 Op 1 dq 
Sw to 4g (— grade + F . (2.8) 
dt dt p oa p dt 


Introducing the last expression into the equation of motion (M’), we finally obtain 


dho __ ds 1 op 
— = T—+ — Pd —q-:(Aq + 4 grad div q). (M’”’) 
dt dt p ot p 
3. The energy equation. From the law of conservation of energy, it can be shown® 
that 
du dp) 1 c 
—— Pp a > + ©. (E) 
dt dt p 
The first term on the left-hand side accounts for the rate of change of internal energy; 
the second term stands for the work required to compress the fluid. The first term on 
the right-hand side represents the heat generated by the viscous forces. The dissipa- 
tion function ¢ is defined by 


d u ov\? Ow \? Ow dv\? 
6 = u[2() +2(=) +2(=) +(= +=) 
Ox oy Oz Oy 02 
Ou Ow? /dv du\? 2... p 
+i — + —}+1—+ —)- — Gr ei. (3.1) 
Oz Ox Ox dy 3 
Finally, the last term on the right-hand side accounts for the heat added to the fluid 


per unit of time per unit of mass. For instance, when all the heat transfer is due to the 
conductivity of the fluid (no external heat sources, no radiation), Q is given by 





QO = p™ div (k grad 7). (3.2) 


The energy equation can be simplified by introducing the entropy from Eq. (2. 3); 
thus we have 
ds 
T—=0 + - (E’) 
dt 
The interpretation of this last equation is particularly simple. The right-hand side 
represents the total heat increase of the fluid, while the left-hand side gives the corresponding 
product of the temperature by the entropy change. By combining the equation of motion 
(M’’’) with the energy equation (Z’), we obtain the following important relation, 


Lho 1 a 
es ate “? 0+ ~¢ pe —q-:(Aq + } grad div q). (E”’) 


dt p ot 


In the literature, the energy equation is frequently given in this last form. 
5 See J. Ackeret, Handbuch der Physik, vol. 7, Berlin 1927, chap. 5, p. 293, or H. Lamb, I.c. pp. 575 
and 637. In the energy equation it is not assumed that yu is constant. 
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Il. VORTEX THEOREMS FOR FRICTIONLESS FLUIDS 


4. The circulation theorems. The circulation is defined by 


r = Pad (4.1) 


where the integration is to be taken along a closed curve formed by fluid particles. It 
is easy to show that the rate of change of the circulation is given by 


aT dq 

eaten pn a Y (4.2) 
dt dt 

Let us combine the last relation with the equation of motion (M). If the viscous 


terms are omitted, we have 


ar 1 dp 
<= - f— srad pa=- § (4.3) 
dt p p 

By means of the identity 


0 = $ dio») = $ pao) + § op (4.4) 


Eq. (4.3) can be transformed into 


aT 
ny - f p(grad p)-dl = - f pate, (4.3’) 


Instead of using Eq. (M), we can use Eq. (M’) and thus obtain 


av 
a = g Trad s)-dl = § Tis. (4.3”) 


Sometimes it is preferable to transform the line integrals into surface integrals with 
the aid of Stokes’ theorem. Thus it follows from the last two equations that 


ar 

= = - $5 [grad p-') X grad p]dA, (4.5) 
rf 

av 

a - $s [(grad T) X grad s]dA. (4.5’) 
C 


We now come to the interpretation of the equations for dI'/dt. When the fluid is 
barotropic, the right-hand side is zero in all of these equations, and the theorem of 
Lord Kelvin is then obtained. The circulation along a closed “fluid line” in a barotropic 
fluid is constant for all time. In particular, when the circulation is zero at a certain 
instant, it will remain so for all subsequent time. By applying Kelvin’s theorem to 
an indefinitely small closed line, Lagrange’s theorem is obtained. 

In the case of non-barotropic fluids, the situation is quite different. The right-hand 
sides in the equations are not zero in general and Kelvin’s theorem does not hold. 
Bjerknes? gave a simple geometrical interpretation of Eq. (4.5). Let us draw equi- 
distant members of the families of surfaces p=const. and p-!=const. and so obtain 
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a series of tubes bounded by these surfaces. The theorem of Bjerknes states that the 
rate of change of circulation per unit of time along a fluid line C is proportional to the 
number of tubes surrounded by C. (In the case of a barotropic fluid the surfaces 
p=const. and p-'=const. are identical.) 

A very similar interpretation can be given to Eq. (4.5’) by considering tubes 
formed by the families of surfaces J =const. and s =const. In the case of a barotropic 
fluid, these two families of surfaces are identical, unless the flow is isentropic in which 
case the surfaces s=const. are no longer defined. Bjerknes’ theorem, in this modified 
form, will be useful in a later part of this paper. 

5. Theorems with respect to the rotation. Helmholtz’s theorem. With the aid of 
Eqs. (2.5) and (C’) it can be easily proved that 

—1 
curl nt = =. — (@V)-q. (5.1) 
dt dt 
Applying the operator curl to both sides of the equation of motion (M) or (M’), we 
obtain in the frictionless case 


d(p~'w) 

aoe id — (oV)-:q = — grad g* xX grad P, (5.2) 
or 

d(p~'w) 

= — (wV)-q = grad T X grad s. (5.2’) 


In the case of barotropic fluids, the right-hand side equals zero. (For two-dimensional 
flows the second term on the left-hand side equals zero because w is everywhere nor- 
mal to q.) Thus, for barotropic fluids, 


d(p~'w) 
dt 





= (p-'wV) -q. (5.2”) 


A geometrical interpretation of the last equation led Helmholtz to the discovery 
of his famous vortex theorems. Vortex lines are material lines. The product of the cross- 
sectional area and of the vorticity w of a vortex filament is constant both in space and time.* 
(Holmholtz unnecessarily restricted his investigations to incompressible fluids.) In 
the case of non-barotropic fluids, (5.2’’) must be replaced by the more general Eq. 
(5.2) and the Helmholtz vortex theorems do not hold any more. Friedman® derived 
certain theorems for non-barotropic fluids which are somewhat analogous to the 
Helmholtz theorems. 

6. The theorem of Crocco and its generalization. In the case of steady, frictionless 
flows the equation of motion (M’’) simplifies to the important relation 


q X w = grad hy — T grad s. (6.1) 


We will see later that for a very important type of flow fo is constant throughout the 
field. In this case Eq. (6.1) reduces to 
* A line which at each point is tangent to the vorticity vector w, at this point, is called a vortex line. 


An infinitely thin tube formed by vortex lines is called a vortex filament. d 
6 A. A. Friedmann, Uber Wirbelbewegung in einer kompressiblem Flissigkeit, Zeitschrift {. Ang. Math. 


und Mech., 4, 102-107 (1924). 
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q X wo = — T grad s. (6.1’) 


This last relation was discovered by Crocco.? When both fo and s are constant the 
right-hand side of Eq. (6.1) is zero and so the motion must be irrotational. The im- 
portance of Eq. (6.1) lies in the fact that it relates the rotation of the fluid to the rates of 
change of ho and s. 


Ill. ADIABATIC, STEADY, FRICTIONLESS FLOWS 


7. General relations. For the flows considered in this chapter the energy equa- 

tions (E’) and (E’’) reduce to 

ca: Ot = (7.2) 

ao ; . ao SC 
where 0/0¢ indicates differentiation along a streamline. Accordingly, both the entropy 
and the stagnation enthalpy are constant along each streamline (but they might vary 
from one streamline to another). Because of its great importance, we shall write out 
the integral of Eq. (7.2) in detail for a perfect gas with constant specific heats. One 
obtains 


1 1 1 Cp Pp ; z 
ho = —PVrth=—_qt+c,T = —q? + — — = const. along a streamline, (7.2) 
2 2 2 ; oe: 
where the equation of state 
p/p = RT (7.3) 


is used. 

The modified Bjerknes theorems simplify somewhat for the flows considered in 
this chapter, because the lines of constant entropy coincide with the streamlines. 
Similarly the generalized Crocco theorem [Eq. (6.1)] simplifies, because the stream- 
lines coincide with both the lines of constant entropy and the lines of constant stagna- 
tion enthalpy. 

An example illustrating these theorems will be useful.* Consider the discharge of 
a perfect gas from a container. We assume that the gas is originally in equilibrium, 
that is, that the pressure po is constant, but do not assume that the temperature 75 
is constant. In order to use Bjerknes’ theorem (in its modified form) we construct 
the net formed by the lines of constant entropy and the lines of constant temperature. 
At the beginning of the experiment the pressure is constant and these lines coincide. 
Thus it follows from Bjerknes’ theorem that dI'/dt=0. However, at a subsequent 
instant, the lines become distinct and so the motion becomes rotational. Let us pro- 
ceed now to determine the rotation. In order to use the generalized Crocco theorem 
we consider only steady state flow (infinite container). According to our energy theo- 
rem, the entropy and the stagnation enthalpy (and consequently the stagnation 
temperature) are constant along each streamline. Furthermore, since fo is a constant, 
it follows from the thermodynamical relation (2.3’) that 


grad hyo = To grad s. (7.4) 
Thus from Eq. (6.1) 


* The author is indebted to Professor H. W. Emmons of Harvard University for this example. 
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T 
axe =(1-—) grad ho (7.5) 


0 
or, after simplifications, 


q X w = 3q*-grad (In 7»). (7.5’) 


We observe again that although the flow originates from a resting gas, the motion is 


rotational in general. 
8. Two-dimensional flow. The continuity equation shows that in this case there 


exists a stream function such that 
u = p'dy/dy, v= — p-'dy/dx. (8.1) 


From the definition of the rotation, it follows that the stream function must satisfy 
the following equation 


O(p Wz) /dx + A(p'y,)/dy = —w (= — 00/dx + du/dy). (8.2) 
In order to determine w we use Eq. (6.1). Because s and fo are constant along a 
streamline, 
Os Oho ; 
qo= T———» (8.3) 
On On 


where 0/0n indicates differentiation normal to a streamline. Both the entropy and 
the stagnation enthalpy are functions of y alone. By using the relation 





: : (8.4) 
= _ eee . 
On . oy 
we find from Eq. (8.3), 
Os Oho : 
w= p r= -—), (8. 3’) 
dy dy 
For a perfect gas this reduces to 
p Os Oho 
wo =———p—: (8.3) 
R dy oy 


Noting that dho/dy and ds/dy are constant along any given streamline, one observes 
that the rotation on each streamline is a linear combination of the density and the pres- 
sure. If ho is constant, throughout the flow the rotation is proportional to the pressure.* 
If s is constant, throughout the flow the rotation is proportional to the density.’ 
(The constant of proportionality is given by the rate of change of entropy or stagna- 
tion enthalpy normal to the streamline.) 

It is of some interest to develop Eq. (7.5’) for the two-dimensional case. Here we 

find that ; 
ee eee: (8.5) 

2 ay 





and the rotation is thus seen to be proportional to pq’ along each streamline. 

Finally we mention that after rather lengthy computations the differential equa- 
tion for y, Eq. (8.2), can be transformed into 

7 K. O. Friedrichs (and R. von Mises), Fluid dynamics, Brown University, Providence, R. I., 1941, 
p. 229. 
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(: “9 a, +(1 “yy _ [om i ( nq q’ =] 6.6 
2? rz a? zy 2 vy — PP ay RR 0 c= . . ) 


When both the entropy and the stagnation enthalpy are constant throughout the 
field, the right-hand side of Eq. (8.6) becomes zero and one obtains the familiar equa- 
tion of a steady isentropic irrotational flow. 

9. Flow around an obstacle with shock-waves. Shock-waves can be included in 
our theory by admitting such discontinuities in the flow pattern as are compatible 
with the laws of conservation of matter, momentum and energy. Thus, the previous 
theory can be applied for flows between shock-waves Fer most purposes one can as- 
sume that the air comes from a homogeneous condition and in particular that ho 
and s are constant far ahead of the obstacle. It follows from Eqs. (7.1) and (7.2) 
that both hy and s are constant at least up to the first shock-wave, and then again along 
each streamline between consecutive shock-waves. Hence w=0 on each streamline up 
to the first shock-wave. In particular, if a streamline is not intersected by a shock-wave, 
w remains zero all along this streamline. From the law of conservation of energy it can 
be deduced that 4» must be continuous across a shock-wave and thus hp must be a 
constant throughout the field. Hence from Eq. (8.3’’) 

o= £ =, (9.1) 
R oy 
and the rotation is proportional to the pressure along each streamline between shock- 
waves. Furthermore it is known that the entropy increases across a shock-wave and 
the increase depends on the magnitude of the shock. Hence 0s/dy is not zero in general 
after a shock-wave and the motion is rotational. Generally speaking there is always a 
sudden increase of the rotation across shock-waves (see Hadamard), and then the rotation 
remains proportional to the pressure (see Crocco). 

10. Flow with axial symmetry. Let the x axis be the axis of symmetry of the flow. 

Then there is a stream function such that 


u = r—'p—dy/dr, v= — r'p"dy/dx, (10.1) 
where 
r=/y?+ 2? (10.2) 


and v is the velocity component normal to the x axis. 
Quite similarly to the two-dimensional case, it follows from Eq. (6.1) that, in the 


present case, 


Os Oho 
w= n(7=-—). (10.3) 
ap ow 
For a perfect gas, we have 
Os oh 
PE i etch (10.3) 
R oy oy 


When fp is a constant throughout the field, one recognizes in Eq. (10.3) a relation 
discovered by Crocco. 
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ON VELOCITY CORRELATIONS AND THE SOLUTIONS OF THE 
EQUATIONS OF TURBULENT FLUCTUATION* 


BY 
3. tea 
National Tsing Hua University, Kunming, China 


1. Introduction. The theory of turbulence, as developed from Reynolds’ point of 
view, is based upon the equations of turbulent fluctuation [1] and has been applied 
to the solutions of various special problems [2, $.4. 5.4; 7|. Owing to present cir- 
cumstances, these papers either have not been submitted to scientific journals for 
publication or are already printed but have failed to appear before the scientific 
public. The theory in its original form and its applications has three apparent diff- 
culties: first, the equations of correlation of the second, third or even higher orders 
constructed out of the equations of turbulent fluctuation contain the unknown terms 
of correlation between the pressure and velocity fluctuations; secondly, there exist 
in these equations the terms of decay of turbulence the values of which have to be 
determined; thirdly, when the differential equations of the velocity correlations of a 
given order are derived from the equations of turbulent fluctuation, the presence of the 
inertia terms causes the appearance of the velocity correlations of the next higher 
order, which are also unknown. This has been pointed out by von Karman and 
Howarth [8] in their theory of homogeneous isotropic turbulence. 

In the present paper we shall show that the pressure fluctuation can be derived 
from the equations of turbulent fluctuation, and is expressible as a function of the 
velocity fluctuation, the mean velocity inside the fluid volume, and the pressure 
fluctuation on the boundary. We shall also show that the decay terms can be put 
into simpler and more familiar forms by kinematic considerations. A general equation 
of vorticity decay will be derived for the determination of Taylor’s scale of the micro- 
turbulence which appears in the decay term; in the case of homogeneous isotropic 
turbulence, this equation was given first by von Karman [8]. To get over the third 
difficulty we shall compare the orders of magnitudes of the different terms in the equa- 
tions of triple correlation. We shall find that the term involving the divergence of 
the quadruple correlation is actually smaller than the correlation between the pres- 
sure gradient and the two components of velocity fluctuation, and can therefore 
be neglected as a first approximation. From this we can also understand why, for the 
flows in channels and pipes in which the mean velocity profile is comparatively steep, 
particularly in the neighborhood of the walls, all the equations of mean motion and 
the equations of double and triple correlation are necessary to describe the phenomena 
of turbulent motions of fluids. On the other hand, as a consequence of the approxima- 
tion based on the fact that the divergence of the quadruple correlation is smaller 
than the correlation between the pressure gradient and the two components of veloc- 
ity fluctuation, we can stop at the equations of triple correlation instead of building 
equations of higher orders. As a matter of fact, for the flows in jets [3] and wakes [4] 


* Received Aug. 21, 1944. 
** Now at California Institute of Technology. 
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where no wall is present, the equations of mean motion and of double correlation are 
sufficient, after some simple approximations to the triple correlations are made, for 
the determination of the mean velocity distribution, and the equations of triple cor- 
relation can be dispensed with. 

From a mathematical point of view the present program indicates that the turbu- 
lence problem can be reduced rigorously to a set of non-linear partial integro-differential 
equations the solutions of which are very difficult to ascertain. In order to facilitate 
the solution of special problems, approximate forms of the integral parts of the equa- 
tions have been developed in a general way. These approximations, however, are only 
valid in regions not too close to the boundary of the moving fluid volume. It may also 
be worthwhile to point out that the unsatisfactory part of the present theory lies 
in the uncertain nature of the correlation integrals, as will be seen presently in §8. 
A better and more accurate representation of these integrals is possible, provided 
more accurate experimental information can be obtained as to the distribution of 
turbulence levels and to the correlation functions between two distinct points in 
general. 

The rigorous way of treating the turbulence problem is probably to solve the 
Reynolds’ equations of mean motion and the equations of turbulent fluctuation simul- 
taneously. This procedure, however, is very difficult owing to the non-linearity of the 
two sets of equations. Hence we have adopted the method of solving the equations 
of turbulent fluctuation by setting up the differential equations satisfied by the veloc- 
ity correlation functions of different orders, a method initiated by von Karman and 
Howarth [8] in treating the problem of homogeneous isotropic turbulence. This 
process of setting up the correlation equations of different orders and seeking their 
solution can be regarded as a method of successive approximation to the solution 
of the turbulence problem; it will be explained in the concluding section of the present 
paper. The correlation functions of higher orders in the various special problems, ob- 
tained by this setting-up process, should be verifiable by direct observation with the 
advance of modern experimental technique; at present experiments have only been 
performed to measure the mean velocity distribution and the second order stress 
tensors in a turbulent flow. It should also be noted that although the equations of 
correlation have a much more complicated mathematical appearance than that of 
the Navier-Stokes’ differential equations from which they are derived, the method of 
Prandtl’s boundary layer approximations can still be used without leading to contra- 
dictions for the particular problems [3, 4, 5] under consideration. 

For the sake of convenience we list below the different equations of motion which 
have been derived heretofore [1]. Reynolds’ equations of mean motion and the equa- 
tion of continuity for an incompressible fluid are given by 

Ou; 





+ UiU;,; = —- he Bait = Tig tvvui, Ui ;=0. (1.1) 
ot p p 

Here, the tensor notation is employed, and U; are the velocity components of the 

mean motion, ¢ is the time, p is the density, # is the mean pressure, v is the coefficient 

of kinematic viscosity, a subscript preceded by a comma denotes the covariant 

derivative, V? denotes the Laplacian operator, and Reynolds’ apparent stress 7’; is 

defined by the relation 
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ri, = — puna, (1.2) 


w; being the velocity components of the turbulent motion. 

The equations of turbulent fluctuation and the equation of continuity for the 
velocity fluctuation w‘, which are the differences of the Navier-Stokes’ equations and 
Reynolds’ equations (1.1), are 
dw, 1 1 ‘ . 
eet Uiw;,;+ w'w;,; + wiU;;= — — 65 — — Tig t+ vV "Wi, ww’; = 0, (1.3) 
ot p p 
where o@ is the pressure fluctuation. From the above set of equations we derive the 
equations of vorticity fluctuation, 


Ph + Uso, + Ui .wi,3 — Ui we, + wlan, + wi .Wi,j — W) We, ; 
t 


1 : 
+ wiQix, ; wi U5 — we Ue = — (81 Te ii) H+ PV wir, (1.4) 
p 


where the mean vorticity 2,;, and the vorticity fluctuation w,, are defined by the equa- 


tions 
Q%. = U © Bees L kyis Wik = Wik — Wk,i- (1.5) 


The equations of double velocity correlation derived from (1.3) are 





1 OT: 1 Part 1 eee 
pow — — (U;, sri, + U4, 574s) — — Urn, + (wis); 
p oF p p 
ae oe v exe tee 
es (@,:We + @,4W:) et Parte V*rik a 2vg™" W;mWk,ny (1.6) 
p p 


where the superimposed bar denotes the mean. The ten equations of triple correlation 
are 





w,w,.w, + U;,;wiw.wit Uy, jwiwyw; + U1, jw wiew, + U' (www); + (wiwiw,.w)),; 





= — — (GO ;wews + ©, WW; + W,,.Wi Ws) 
p 


1 ue aes, 
+ = (745, j7Kt He 72K, Tu HE TA, Tix) H Vg" (WiWEW2) ,mn 
p? 





— 2vg™"(Wi,mWknWi + Wk,mWinWi + Wi,mWinWk). (1.7) 


2. The pressure fluctuation. Let us take the divergence of the equations of turbu- 
lent fluctuation (1.3). Because of the equation of continuity satisfied by w‘, the pres- 
sure fluctuation @ satisfies the following Poisson’s equation: 


1 ee cncentne 
— Va = — 2U" wm + (Ww Ww" — ww") mn: (2.1) 
p 


Since any two successive covariant differentiations are commutative in a Euclidean 
space, the gradient of the pressure fluctuation ®,, also satisfies a Poisson’s equation, 
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1 es 
— Va. = — 2(0U" , Wm) k + (W"™w" — WwW") mnk (2.2) 
p 


The general solution of (2.2) can be written in the form, 


1 1 
—.= fff (U'* a" «)’ 2 —av’ — SIS (w’ ™w’™ — w'™w!")’ snk —av’ 
p 2 
1 do’, 0 1 
1 ff jt ay Sy” 
we r on’ on’ \ r 


where the integrations extend over the whole region of the moving fluid, the first 
two integrals represent the particular integrals, and the third represents the comple- 
mentary solution which is a harmonic function expressed in terms of the boundary 
values of itself and its normal derivative; x’' are the coordinates of a point P’ which 
ranges over the region of the moving fluid, 7 is the distance from P’ to the point P 
with coordinates x‘, dV’ is a volume element, dS’ is a surface element, 0/0’ denotes 
the normal derivative, and the primes on the various quantities on the right side of 
(2.3) indicate that these quantities are to be evaluated at P’. We shall see finally 
that the surface integral in (2.3) can be, neglected for points P where o,, is defined 
and which are not too close to the boundary of the moving fluid. 

We now let both x‘ and x’‘ represent rectangular cartesian coordinates, and let & 
denote the difference vector of x’‘ and x‘, i.e., 


tim gi — xi, (2.4) 








Covariant differentiation then reduces to ordinary differentiation, and the difference 
between covariant and contravariant tensor character disappears. Hence £& is equal 
to £; and the distance r between P and P’ is given by 


2 = EiE%. (2.5) 


The element of volume dV’ is equal to dé'd&*d&. 

The solution (2.3) clearly shows that besides the harmonic function expressed as 
a surface integral on the boundary, the pressure fluctuation at a point P, and its 
gradient, are determined by the turbulent velocity fluctuation w‘ not only at P but 
also everywhere within the fluid. However, due to the factor 1/r in the integrands 
the effect of the velocity fluctuation at distant points P’ on the pressure fluctuation 
at P gradually dies away as P’ recedes farther and farther from P. 

3. Velocity correlation between two distinct points. The partial differentiations in 
the integrand functions in (2.3) are taken with respect to the coordinates x’* which 
are independent of x‘. Hence, if we multiply (2.3) by the velocity fluctuation w; at 
the point P, we obtain the correlation between w; and @, at the same point P: 


1 Sees 1 1 ialeaitaibaaiete 1 
—@ .W; = —ff [U’™ ,(w’"™ wi)! om} _ dV’ + he (w’™w'"w;)’ mak — dV" 
p 2x r 4 r 
1 
—ff {- = — “(-) dS’. (3.1) 
“we r On’ 


We shall neglect, however, the surface integral in the above equation on the 
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ground that the correlation g’,w; is small provided that the point P where the correla- 

tion @,w; is under consideration is situated not too close to the boundary. This 

condition limits the present theory to regions where free turbulence predominates. 
Likewise, under the same condition of approximation the correlation function 


pa ,wyw; is given by 


5 eee 1 ee 1 
—@,.wWiw; = —{f [U’™ ,(w’™ ww)’ m |’ x — dy’ 
p 2 r 


1 —_ Ls See 1 
- —fff [w?™w’™ ww; — w ™ www; |’ mak — AV’. (3.2) 
An r 


If we solve for o from (2.1) and form its correlation with w;, we find, to the same 
order of approximation, that 


cof s eee. 1 1 
— ns fff L ; m n(w’"w;)’ m At dV’ + pack fff (w’ may! ™w;)! mn sg dV’. (3. 3) 
2x r Ar r 


p 





In the three equations (3.1), (3.2) and (3.3) we recognize three types of functions, 
namely, w'™w,;, www; and w*ww;, and w'™w'"w,w;; they are, according to Taylor 
[10] and von K4rmé4n [8], the velocity correlations between two distinct points P 
and P’ of the second, third and fourth orders respectively. They are usually functions 
of both the coordinates x‘.and x’* and probably also of the time ¢. The double correla- 
tion function w’"w; between P and P’ has been measured extensively for isotropic 
turbulence by several authors [11, 12]; for flow in a channel [13] and in a pipe [14], 
they have been recorded only in a number of isolated cases and only within limits. 

It has been observed that for isotropic turbulence w’*w; vanishes very rapidly 
for large values of the quantities £* defined in (2.4). This must also hold true for the 
other two correlation functions w’™w’*w; and w'*ww;, and also for other types of 
flow ; furthermore their derivatives with respect to £* should all approach zero rapidly 
with increasing £*. 

On the other hand, the quadruple correlation w/"w’"w,w; between the points P 
and P’ does not necessarily vanish when P and P’ are widely separated, for the aver- 
age values of both w’"w’" and wyw; over a period of time 7 are themselves not sepa- 
rately equal to zero in general. Hence as an analogy to the velocity vector W‘, we 
may separate the product w,w; into two parts, the correlation ww; and a symmetric 
tensor “;; the time average of which vanishes, * 

















WwW; = WwW; + UH, (3.4) 
ee 1 rtr/2 
i uyidt = (). 
i r—r/2 


An analogous relation holds good for w™w’. The quadruple correlation w/™w/*wjw; 
between the points P and P’ consequently becomes 





w ™w' "ww; = w ™y'™ W1W; + ul 4 1;. (3.5) 
As P and P’ recede farther and farther from each other, the correlation function 


* The author wishes to express his gratitude to Mr. S. L. Chang for pointing out relation (3.4). 
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u’™"u,;, which behaves like w’"w;, will tend toward zero as a limit. Substitution of 
(3.5) into (3.2) yields 


ee oe 1 seiaiecoitca 1 
— BWW; = —{ff [U'™ ,(w’™wywi)’ m |’ — dy’ 
p 2r r 


+2 ify (Ta) ua — BV" (3.6) 


If we substitute in the equations of the double and triple correlations (1.6) and (1.7) 
for p-'w ,w,; from (3.1) and pe ,w,w; from (3.6) above, we obtain a set of integro- 
differential equations for the mean velocity, the double and triple velocity correlations 
of a turbulent flow at a point P being the dependent variables with the velocity correla- 
tions between two distinct points w’"w; and w’*w,w,; as kernels. This set of integro- 
differential equations is too complicated for solving special problems, so we shall 
presently develop approximate forms of the integral parts of the equations in a gen- 
eral way. 

It should be noted that for homogeneous isotropic turbulence the following rela- 
tion between the triple correlations holds [8]: 











W' mW ,Wi = — W pWmWn. (3.7) 


4. Conservation relations satisfied by the velocity correlations. The velocity fluc- 
tuation w’" at the point P’ satisfies the equation of continuity w’,, =0. Let us multiply 
this equation by w; and average over an interval of time r. Since P and P’ are inde- 
pendent, we obtain the conservation equation for the double correlation w’"w; be- 


tween P and P’, 





(w’"wi)2 = 0, (4.1) 
Ox’" 
where the coordinates are still rectangular cartesian, and the subscript x indicates 
that the variables x* are to be held constant while the differentiation is carried out. 
Instead of x‘ and x’‘, we can use the coordinates x‘ and £', i.e., we transform from 
the old variables x* and x’‘ to the new variables x‘ and £* by means of the equations 





ai = gi, xt = git FF, (4.2) 
In terms of the new coordinates x‘ and &*, Eq. (4.1) becomes 
(ww). = — (Ww). = 0. (4.3) 
0x" og" 


For the sake of simplicity we shall drop the subscript x in (4.3); it will be understood 
that the variables x* are regarded as constants during the differentiation. Hence we 
can write the divergence equation (4.3) in the covariant form, 


(w’"w;).» = 0. (4.4) 


Similarly, from the equation of continuity for w‘ at the point P, we have 


te) 
ar (w’,Ww*) s = 0. 
Ox* 
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In terms of the new coordinates x‘ and £‘, this relation becomes 


fe) —— fe) —— 
Po (w’,w*): — ~ (w’,w*): = 0. 

x ‘ 
After changing the variables from x‘, x‘ to x‘, £' it can be seen that w’"w;, considered 
as a function of x‘ and £‘ rather than of x‘ and x’, varies slowly with x‘ but rapidly 
with £' for points not too close to the boundary of the fluid volume. Hence, as a first 


approximation the equation of conservation for the double correlation between P 
7 , iy sage 
and P’ in the index 7 is given by 


(w’,w'),; = 0. (4.5) 


We note that to the first approximation the correlation function w;/ w; satisfies the 
conservation equation symmetrically with respect to the indices i and k. 

Likewise, the other two correlation functions w’™w'"w; and w'"w,w, between P 
and P’ can be shown to satisfy the following relations: 


(w'™w'"wi),; = 0, — (w/™wiwe) n = 0. (4.6) 


The first equation in (4.6) is derived by an approximation as was (4.5); the second one 
is rigorous. We must not forget that all the covariant derivatives in (4.5) and (4.6) 
are taken with respect to the variables £‘, the coordinates x‘ being held constant. 

It is obvious that since the coordinates x‘ of the point P are regarded as constants 
under the integrations in (3.1), (3.2) and (3.3), the covariant derivatives with respect 
to x’* in the integrand functions can all be replaced rigorously by covariant deriva- 
tives with respect to the variables &*, because of the equations of coordinate trans- 
formation (4.2). For example, (3.1) then becomes 


_ 1 aa 1 
—@,W; = fff [U’™ ,(w'" wi) mx — dV’ 
p Qn r 
+ SS for Wi) smn — av" (4.7) 


The other two integrals (3.3) and (3.6) can be altered analogously 

5. Correlation integrals between the pressure gradient and velocity fluctuations. 
Let us examine the integral (4.7) more closely. In the integrand function of the first 
integral on the right hand side, U’",, is a more slowly varying function of £* than 


its factor w’"w;, both functions being regarded as functions of x* and &*. Hence, we 
expand U’”,,, at the point P’ in a multiple power series in £', 


1 etiyy u 
S 








au’™ = au™ 


sree = — to 


Ox” Ox” pee 


~ ghgle... gle, (54) 


Ox dx? ~~ + Axsdx” 





Substitution of (5.1) into (4.7) would yield a series of integrals which would be too 
complicated for any practical application. But if we neglect the higher order terms, 
in (5.1), then we have as a first approximation to (4.7), 

1 


— (0; {Wk “+ O..w y;) = aia a" nia ™ ~T bik, (5.2) 
p 
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where the functions a",;, and 6, are defined by 


1 a ae a 1 
a" nik = —fff [(w’™ 2;) mk + (w’" w,) mi| — dV’, 
2r r 
1 ee aid pet oosie t 1 
bs a tone fff [(w’ ™w’"w;)  mnk + (w’ ™w’" wy) mni| —_i. (5.3) 
4a r 


Owing to the conservation relations (4.5) and (4.6), the above two sets of functions 
also satisfy the following divergence conditions: 


a" nik - 0, g**a" nik - 0, gb = 0. (5.4) 


Of these three conservation relations, the first follows from the rigorous continu- 
ity equation (4.4) and is hence exact, while the other two follow from (4.5) and the 
first equation of (4.6) and are hence approximations. The nature of the functions @"nix 
and 6;; will be discussed in §8 below. 

Because of (5.4), contraction of (5.2) by means of g** yields, 

ee ea et 

—o,w' = — (ow); = 0. (5.5) 

p p 
This result is consistent with the correlation (3.3). For we may substitute the series 
in (5.1) into (3.3) and preserve the largest term; but the latter is smaller than the 
first term on the right-hand side of (5.2) by a factor of \ which is Taylor’s scale of 
micro-turbulence [10, 8]; the second term on the right-hand side of (3.3) is also smaller 
than };, by an analogous factor. Hence the approximate form of (3.3), to the same 
degree of accuracy as in (5.2), is 





1 
— aw; = 0. (5.6) 
p 
This relation has also been proved to hold true for isotropic turbulence by von 
K4rman and Howarth [8]. 

By a similar process, we find from (3.2) that the triple correlation between the 
pressure gradient and two components of the velocity fluctuations is, to the same de- 
gree of approximation, 

1 ry 
— (@ ;W,W, + @ , WW; + @® WW) = b" mikil ™ in + Cikly (S. 7) 
p 





where the forms of the tensors b", jx: and ¢jx: are given by respectively by 


1 i acta Pa tat: eat 1 
Bb" miki = >, fff [(w’" wwe) me + (w’"w,W1) .mi + (w’" ww) me | « dV’, 
LT 


1 5h eanhh heme ot ws 1 
second 4 Sf [(10’™" 105%) ,mnt + (ul ™" UK), mni + (u’ ™41;) mnk | Pal (5.8) 
T 


Because of (4.6), the functions b",,;.; satisfy the rigorous conservation relation, 
b" nixt = O. (5.9) 


We shall discuss the general behaviour of the functions b"mie: and cix: in §8. 
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6. Terms involving the decay of turbulence in the equations of double and triple 
correlation. To determine the terms in the decay of turbulence, it is necessary to 
know explicitly the double and triple velocity correlations between two adjacent 
points. Physically, the correlation functions between two near-by points must satisfy 
two conditions: first, they should become the velocity correlations at one point when 
the two points coincide; secondly, they should degenerate into the isotropic corre- 
lations when the flow obeys the condition of isotropy. By two adjacent points we mean 
that expansions of the double and triple correlation functions in terms of the coordi- 
nates £* stop after the second and third powers of £*, respectively. Furthermore, since 
only approximate expressions of the decay terms are required, conservation equations 
in the forms (4.4), (4.5) and (4.6) will suffice for the present purpose. In view of the 
property that the double correlation w,w’, satisfies the conservation relations (4.4) 
and (4.5) symmetrically with respect to the two indices 7 and as a first approxima- 
tion, it should also satisfy the supplementary condition that its expansion be sym- 
metrical in the coordinates of P and P’. 

The second order velocity correlation between two adjacent points that satisfies 
the above two conditions and the supplementary condition of symmetry can be ex- 


panded into powers of £ in the form, 





shinies A Ouk ee 
oor, me agt¢—— £.t, By mén = ___. EMER 
wid, = 49 {— bibs + > Bank ™§* + Ru(1 +e ) 


G 1 
- ne (Riké'Ex + Reik'ts) + tne Exe jpmntit'E"E" + - , , (6.1) 


where g is the mean magnitude of the velocity fluctuation, or the root-mean-square 
of the velocity fluctuation, defined by 








gq’? = ww’, (6.2) 
and Rj, stands for 
3 
Ru = — WiWr. (6.3) 
q? 


The function \ is Taylor’s scale of micro-turbulence, both g and A being functions of 
the coordinates x‘ of P; A, Bun, Cun, G and Ejxjimn are all independent of £*. The co- 
efficients B,,, and Cy», are symmetric in m and n; Ejkj1mn is both symmetric in 7 and k 
and in the last four indices j, 1, m and n, but is not symmetric in any one index of the 
first set of two and any one in the last set of four, e.g., it is not symmetric in 7 and j. 
Hence this tensor has 6X15 =90 independent components. 

The form given in (6.1) for the correlation tensor wav’, between two adjacent 
points is the most general linear combination of the products of the tensors ££, dix 
and w,w;,. The functions A, By, C;, and G will be assumed to be constants; it is not 
necessary to know the exact nature of the separate components of £;;jimn for our present 
purpose, but we shall assume for the time being that the invariant E =g‘*g?'g™"E jx jimn 
is constant. 

A question naturally arises as to whether the functions g? and 2? in (6.1), which 
vary with the coordinates, should be replaced by expressions which are symmetrical 
in the coordinates of P and P’. However, this is not essential, for both g? and \? vary 
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much more slowly than w,w’, as a function of £*; since P and P’ are close to each other, 
we may use their values at P as an approximation. Nevertheless, one must be careful 
with this approximation, whenever differentiation with respect to x’‘ is involved. 
The function w,wv’, must satisfy the equation of continuity (4.4), or (4.5). By 
setting the coefficients of —* and £4¢'" equal to zero separately, we find that 


2Abix + Biz + R'C ed SGRix — 366%. = 0, (6.4) 
g**Eizeimn = 0. (6.5) 


Since Eq. (6.4) is symmetric in the indices i and k, and since C,. has been assumed 


to be a constant, we must have 
Gane —_— CSan (6.6) 


where C is a constant. On the other hand if Ca, depends upon the correlation tensor 
ww», then it is possible to have the more general solution Can = —C8mn+DSmn, where 
Smn is the inverse matrix of Rn» defined by S';R*,=6*;. For the sake of simplicity, we 
choose D to be zero for the time being. Obviously, the number of independent equa- 
tions in (6.5) is 30. 

In order to give a simpler appearance to the final forms of the decay term in the 
equations of double correlation and of the equation of vorticity decay, we put 


A=1+4G, C = $(k — 46). (6.7) 


The first equation amounts to a change of the factor A, this factor being arbitrary; 
the change makes \ assume the same numerical value as Taylor’s scale of micro- 
turbulence, when the correlation tensor obeys the condition of isotropy. The second 
equation in (6.7) only defines C in terms of a new constant k. Utilizing relations (6.4), 
(6.6) and (6.7), we put (6.1) into the form 


2 


WW = Wwe + = faa + 4G)EE — 4[(2 + SG)r? — 4(k + 11G) René "€" J5ix 


— (k — 4G)? Riu — G(Rik'& + Riié'Es) + Exnjamnk§'E"E" + + - ; » (6.8) 


1 
“4? 
where the tensor Ejxjmn Satisfies the thirty linear equations (6.5). We shall see pres- 
ently that, with the form of ww’, given in (6.8), only the constant k will appear in 
the term involving the decay of turbulence (6.14), while only G will be present in the 
equation for the decay of vorticity (7.11). 

For isotropic turbulence we have ww, = }q75;x, and it is easy to verify that in 
(6.8) the terms in £;£, and r? coincide with terms in the isotropic correlation tensor 
according to von K4érmaén and Howarth [8]. The validity of formula (6.8) and its 
properties can be subjected to experimental verification. 

For the triple velocity correlation www’, between two neighboring points, we 
have to assume a form which degenerates into wjwjv, when the points coincide and 
becomes the triple correlation for isotropic turbulence when the condition of isotropy 
is satisfied by the flow. Since the expansion of the triple isotropic correlation function 
begins with the third powers of £‘, as shown by von K4rman and Howarth [8], the 
same must hold for the present general case. This expansion must satisfy the equation 
of continuity (4.6), and the final result obtained is 
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Fq 
313+/3r3 





~ Jie 
WWW, = WWW, + 


5 
[Et jE, = 5 inks + 6 éi)r? + 5; ;é4r?| aa SP? « (6.9) 


This equation tells us that up to this degreé of accuracy the correlation function 


w,w,w', is the sum of w,w,w, and an isotropic correlation tensor; similarly, we have 
7 iWjwk ’ y ’ 


ae Fq 
Wi Wm W' n = WiWmWn — - = 
3!3+/3r8 














5 
[2E:Emén mod 5 lomibn + Snitm)r? + Smntir? | + a sli (6. 10) 


In the above expression the relation (3.7) for isotropic turbulence has been utilized. 
In the expansions of (6.9) and (6.10), we have introduced the further assumption 
that the triple correlation h can be expressed by [8] 
F 
r 
3!X3 





3 (6.11) 


k= : 
where X is Taylor’s scale of micro-turbulence and F is a numerical constant which 
may be different for different flows. This emphasizes the point that this length A plays 
an important role, not only for double but also for triple correlations as well. The 
validity of this point should be tested experimentally. 
Differentiating the correlation function (6.8) with respect to x”, we obtain 
a Ow’, een 
WiWe = Ww = —— (wiw's) = 


ax'* ‘Ox'* age 








9 P 


£. ue + 4G) (5iske + Secki) — [(2 + SG)E. — 3(R + 11G) René” 5x 


— 4(k — 4G)E Rix — G(Riske + Ritdesk! + Ricki + Rewdisé') 
1 
+ ary Ex rjimstit =" eye \ ’ (6.12) 
3!? 
and furthermore, under the same approximation as in (4.6) where 0( );/dx‘ is neg- 


lected, we get 


Ow; OW: af dw 3? intel 
-= | —I Ww; - a. pao Tye & 
Ox? Ox* Ox? Ox */ Jz—0 OEP OES t<0 
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re {$(1 + 4G) (iden + Suedip) — [(2 + SG)b.p — (RB + 11G)Rey ]Oin 














— 1(k — 4G)b.pRix — G(Risdep + Ripdts + Resdip + Ripdis) }- (6.13) 


Hence the term that represents the decay of turbulence in the equations of double 


correlation (1.6) is equal to 





dw; Owe 2v (k — 5)a°en + 2vk (6.14) 
— = — —(k — 5)q°gix + — wim. 
Ox™ Ox” 3d? . ? ; 





2vg mn 


If we differentiate the triple correlation (6.9) with respect to the coordinates x’, 


the result is 
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0 — =- Ow’; 0 ——,- 
Pr WiW.W 1 = WiW; PPT y = aem (ww. ») z 
Similarly, to the same order of approximation as in (6.13) the following relation is 


true: 





F) ow) Ow, Ow’, Ow; Ow) a 














(w; ww’ 1) 2 


—— WW - 


ee eae PER, Wk 
Ox” dx’™ Ox" dx'™ dz* dz'* og"0g" 


This equation and formula (6.9) then yield 














” 
Wi; 








OW; Ow, Ow; Ow, 3? aioe 
fonts dp =— WWW 1 = 0. (6.15) 
ax" ax™ ax"ax” aEmag ad 


Cyclic permutation of the indices i, k, ] in (6.15) gives rise to two similar relations; 
the sum of the three is identically zero, which shows that the term analogous to the 
decay of turbulence in the equations of triple correlation vanishes in general: 








Qvg™" | Wim Wk,n Wi + We, mWinWi + Wi,mWinWe| = 0. (6.16) 


7. The equation of vorticity decay. Since Taylor's scale of micro-turbulence \ 
plays a very important role in the decay of turbulence, it is necessary to find the equa- 
tion which governs the behaviour of this fundamental length. This equation is pro- 
vided by the decay of vorticity. The root-mean-square of the vorticity fluctuation 
(w?)/2 satisfies the equation 


= Lymrgntonwnn (7.1) 


where &m, is the antisymmetrical tensor defined by (1.5). It is not difficult to derive 
the equation satisfied by w from (1.4) directly. However, this procedure would be 
too lengthy and we shall pursue an alternative course. 

We notice that 


1 t+r/2 
gg" a mnW pe = —f (Wm,n -_ Wn,m)(w™ .g"* oe w" pg™?)dt 
T t—r/2 


= 2(WmnW™ 2f"* — W" mW n). (7.2) 


On the other hand, to the same order of approximation as in (6.12) and (6.13), the 
following expressions are true: 








saudi seiandoneli aici 0? 
BW, nW™ 2 = — (V7Wmw’™)em0, slit alta tna ( wna") = 0, 
dg™dE" i=0 
where V*; stands for the Laplacian operator in the variables £*. It then follows that 
w? = — (V*pwmw’™) eno. (7.3) 

Our next step is to derive the differential equation satisfied by (V?;w,,w’™)¢=0. 

From the equation of turbulent fluctuation at the point P’, which can be written 
in the form 

Ow’: 1 


' Jere 
“-- + U'iw' ¢.; + w iw’ p,; + wiU':.; i: a: ot — <=} ray; + vV'2w'., (7.4) 
Ot Pp p 
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we derive the equation satisfied by the general double correlation function: 


DO eneye — eee ——7-* — 
— ww’, + Ui(wyw',),; + Ul i(wyw's)’ 5 + (wi wiw's),5 + (w'wiw's)’,; 
+4. w ,wiu;,; + ww iU" g.; 
a woe at se : 
= — — (ww’;).s — — (@w;)’ ct + V?(wiw's) + V'*(wiw's), (7.5) 

where the covariant derivatives ( ),;and ( )’,;are taken with respect to the variables 
xi and x’), respectively. 

In Eq. (7.5) we next replace x‘ and x’‘ by the two new sets of variables x* and & 
by use of (4.2), and neglect terms involving the partial derivatives with respect to x‘ 
when £ are held constant, except for the term U‘%0(_ );/dx‘; this exception is made 


because U’ is large when compared with w*. Since we are only interested in the cor- 
relation functions for two adjacent points, we can write 


a at eee se 
i aes Jony’ - 
(w’iw,;w',)’,; = 3E (wiw’ :we) 2, 


as in the case of isotropic turbulence (3.7). With all these approximations in view, 
Eq. (7.5) in rectangular coordinates then becomes 


O Vek. AO eam he et © WE vi eoacaed 
— ww’, + Ui — (wyw's), — Ui — (wiw's) 2 + U"' — (wyw's)s 
ot Ox! Of? Og? 











x 
a a * - OU’, 
— — (www, + wiw' wry) 2 + ww? teh... Adieorg 
qe ? Of? 
~ — (oa) ~ an Teale + 270% ‘k) (7.6) 
2 ——(ee ie er ee VV": WiW &k) z. bas 
p og p og* 


For two adjacent points the power series expansion of aw’, in £* is in odd powers of £'. 
Hence, by interchanging the two points P and P’, we should have 


ow; = — ww. ( 





~I 
~~ 


“I 


Consequently (7.6) is essentially symmetric in the indices i and k. 
Next, let us contract the indices i and k in (7.6). As in (4.4), ww’, should satisfy 


rigorously the equation of continuity, 
° @w), =0 (7.8) 
— (ow’*), = 0. ; 
o¢* 


The result of this contraction then becomes 





a Py phlei el naa F) 
— ww'* + Ui — (ww'*); — Ui — (uw'*), + VU" (w,w’*), 
Ox? 0&? 0&? 


ot 

ae er i “eo 
—2 oe (wiw,w’*), + ww! —— + ww’? —— = 2v9?;(urw’*) z. (7.9) 
7 
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Let us operate upon (7.9) with the Laplacian operator Y*; and then set £*=0, 
denoting ( ):~0 by ( )o for smiplicity. The resulting equation is, 
0 o? 


pow sae 0 ee 
— (V*,w.w’*)o + Ui — (Vw w’* 2U? mg” ( ww) 
at gWEW “)oO aah EWk )o + 4 agiag" k : 








a ee alin ial 
—2 (v' aE wai) + 2U x, (V7? ww’ *)o = 2v(V4¢w,w' *)o. (7.10) 
0 


It is to be noted that in the above equation we have neglected the term w,w'V?,U*, ;, 
which is smaller than the term U;,;(V?¢ww’*)o by a factor which is the square of the 
ratio of \ to a macroscopic length. Equation (7.10) also follows from the equation of 
vorticity fluctuation (1.4) directly as mentioned before. 

By substituting into (7.10) the explicit forms of the correlation functions ww’; 
and w,w,w', for two adjacent points given in (6.8) and (6.9), respectively, and then 
setting £'=0, we obtain the equation of vorticity decay, 





_ 8/9 8 (¢? 14G —— 86 2 _ ¢ 
5— <) + sui — C)-—vuwer- — fe Sk, (7.11) 
dt \ ? Oxi \ 2 d? 33 3 3 4 
in which £ is defined as before, 
E = gitgigm” Et jimn- (7.12) 


We assume that both E and F are constants which may be different for flows with 
different Reynolds numbers. In deriving equation (7.11), the equation of continuity 
U’,;=0 for the mean motion has been utilized. It is also readily verifiable that (7.11) 
agrees with von Karm4n’s equation of vorticity decay for isotropic turbulence [8]. 
8. Nature of the correlation integrals and the final forms of the dynamical equa- 
tions of correlation. Up to the present the only remaining uncertain quantities in the 
equations of the double and triple correlations (1.6) and (1.7) are the correlation 
integrals, a@"mix, and by, in (5.3), b mie: and cya: Of (5.8), and the quadruple velocity 
correlation w,w,w,w,. Let us’ examine the correlation integrals first. The function 
@"mik defined in (5.3), for example, would be uniquely determined if the double cor- 
relation ww’, were known. But unfortunately the equation of continuity (4.4) and the 
general dynamical equation of double correlation (7.6) are insufficient to yield a defi- 
nite solution for w,w’,, because of the presence of the triple correlation w,w,w’, in 
(7.6). 
On the other hand, although the integrand functions of the four kinds of correla- 
tion integrals are not known, we are dealing primarily with the integrals themselves 
and they can only vary slowly with the coordinates involved. This argument can be 
understood, if we recall that the correlation functions w,w’,, wa,w',, ww’ ,w™ and 
U,U'm, under the integral signs only change slowly when both the point P and the 
point of integration P’ undergo a rigid body translation, and that they vary rapidly 
when the relative displacement of the two points changes. This rapidly varying part 
of the functions is integrated away, leaving the slowly varying part behind. The 
neglecting of the term 0(w’,w*)¢/dx‘ against 0(w’,w*)./d£' in (4.5) also follows from 


this interpretation. 
There is another mathematical reason for the fact that the four kinds of integrals 
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are slowly varying functions of the coordinates. If, for instance, we differentiate with 
respect to x* the quantities c;,,; defined in (5.8), we find that 


OCjK1 1 j fe) as —— a 1 
‘ae ae Sf a [( u’ as 53) i + (u! 1441) mni + (u!™™147;)  mnk | we, 
Ox* 4 lax r/i 


re) eee ae ee 1 
— — ( [san + (1! ™"144-1) ni + (u’™"141;) mnk | ~) av (8.1) 
Ss z 
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The first part of the integrand function is small when compared with the second, and 
the second can be transformed into a surface integral on the boundary of the fluid by 
means of the usual divergence theorem of vector analysis. If the point P is not very 
close to the surface, this surface integral is negligible on the ground that the correla- 
tion function u’™"4,;, and its derivatives between P, the point in the interior of the 
fluid, and P’, the point of integration on the boundary, are negligible. 

Since the correlation integrals are slowly varying functions of the coordinates, we 
shall expand them as powers of the coordinates used in the special problems to be 
solved. From kinematic considerations, the integrands of the integrals may further- 
more contain powers of g, the root-mean-square of the velocity fluctuation, as fac- 
tors. Both theory and experiment at present do not assure us of the exact dependence 
of this factor. Nevertheless, so far as the mean velocity distribution is concerned, this 
uncertainty is probably not important, as we shall see in the problem of pressure flow 
between two parallel infinite planes [9]. 

By substituting into Eqs. (1.6) and (1.7) the approximate forms of the four cor- 
relation integrals from (5.2) and (5.7), and the decay terms (6.14) and (6.16), we ob- 


tain finally 
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aa a pe pnd Men ae = 
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Wee 


= BU", — Cnr + = (74, sree + Tie, gris H Tit, jTKs) + vg" (WiWeW:) mn. (8.3) 
0 

In the second set of equations we notice that the term involving the quadruple cor- 
relation is actually smaller than the terms b",:.,U™,, and ¢;x: which form the correla- 
tion between the pressure gradient and two components of velocity fluctuation. This 
is due to the fact that the term (wiw,w,w,),; is equal to a velocity fluctuation raised 
to the fourth power and divided by a macroscopic length, while on the other hand ¢;x: 
is, from its definition (5.8), of the order of a velocity fluctuation raised to the fourth 
power and divided by a length which has the same order of magnitude as Taylor's 
scale of micro-turbulence. The permissibility of neglecting the terms (w‘w,w,w,),; and 
p~*r?; ;7x: as a first approximation, for instance in the problem of pressure flow be- 
tween two parallel infinite planes [9], can be regarded as a justification of the above 
approximation and its associated interpretation. 
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We must not forget that the other dynamical equations necessary for the solution 
of a turbulence problem are the equations of mean motion (1.1) and the equation of 
vorticity decay (7.11). 

9. Conclusion and summary. It is now not difficult to see that the foregoing de- 
velopment is essentially a method of successive approximation to the solution of the 
turbulence problem. In the initial approximation we have the well-known Reynolds’ 
equations of mean motion which contain the unknown apparent stress. From the 
mathematical point of view the momentum and vorticity transport theories connect 
this stress with the mean velocity by physical arguments, in order to make the mean 
velocity distribution determinate. 

The next approximation in solving the given turbulence problem is to use the 
equations of mean motion and of double correlation by making certain approxima- 
tions to the triple velocity correlation in the equations. This procedure has been fol- 
lowed in the determination of the velocity distributions in jets [3] and wakes [4], 
where free turbulence predominates; for the triple correlations we use their values 
at the centers of the flows as an approximation. The mean velocity distributions thus 
obtained agree with the experimental observations very well over large portions of 
the flows. 

In the third approximation to the solution of the problem we have to solve the 
equations of mean motion and of both the double and triple correlations simultane- 
ously by assuming approximations for the quadruple correlations. It is obvious that 
this process of forming the differential equations of the correlations out of the equa- 
tions of turbulent fluctuation can be generalized to higher orders. Fortunately, as in 
the problem of pressure flow through a channel [9] where a wall is present, we can 
stop at the equations of triple correlation and neglect the quadruple correlations as an 
approximation, so that the solution of the problem is not too unnecessarily compli- 
cated from the theoretical point of view. As we shall see, the solution of this particular 
problem holds true in all parts of the channel, if all the equations of mean motion and 
of double and triple correlation are used. On the other hand, the solution for the mean 
velocity based upon the equations of mean motion and of double correlation by using 
the value of the triple correlation in the center of the channel as in jets and wakes, is 
only valid in the central part of the channel, and fails when the wall of the channel is 
approached. This brings up incidentally the important role played by the triple cor- 
relation in such problems. 

In order to see more clearly how the equations of double and triple correlation 
in the forms (8.2) and (8.3) and the equation of vorticity decay (7.11) are derived 
from the equations of turbulent fluctuation, it might be of interest to sum up the 
conditions and approximations under which they are valid. They are listed below: 

(1) The velocity correlation between a point in the interior of the fluid and another 
on the boundary is negligible. This excludes the immediate neighborhood of the 
boundary of the fluid as a region of application of the theory. 

(2) The variation of the mean velocity is small as compared with the correlation 
function between two distinct points when the relative displacement between the 
points changes, so that the higher order terms in the series (5.1) and similar series 
may be dropped. 

(3) The second and third order velocity correlations between two adjacent points 
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are expansible as power series in £*/A with the terms that do not contain £ propor- 
tional to the Reynolds stress at the points. This brings out the point that Taylor’s 
scale of micro-turbulence \ plays an equally important role for both the double and 
triple velocity correlations. 

(4) The slowly varying nature of the functions a"niz, Dix, b%mixr and Ci. with the 
coordinates, and its physical interpretation, have been explained in the preceding 
section. 

With the advance of modern experimental technique the above four conditions 
and their theoretical consequences, as presented here, can all be tested by direct ex- 
perimental observation. The less certain part of the theory lies probably in the dis- 
cussions in §8 of the slowly varying nature of the correlation integrals with the 
coordinates; this perhaps could be improved if more accurate experimental evidence 
were available. 

This paper was written in China before the author’s arrival in this country in 
November 1943, and has been supplemented and revised in Pasadena. It is a great 
pleasure to the author to thank Dr. R. A. Millikan and Dr. Th. von Kérman for the 
opportunity given him to work at California Institute of Technology. He is also grate- 
ful to Dr. von K4rma4n for his interest in the problem and for many helpful criticisms 
and discussions. To Dr. C. C. Lin, who collaborated with the author during the early 
days of the development of the theory some years back in China, the author wishes to 
express his indebtedness for many discussions of the present paper. 


REFERENCES 


[1] P. Y. Cou, Chin. Jour. of Phys., 4 (1940), 1-33. 

[2] C. C. Lin, Pressure flow of a turbulent fluid through a circular pipe (unpublished). 

[3] C. C. Lin, Velocity and temperature distributions in turbulent jets, Sci. Rep. of National Tsing 
Hua University (1941) (already printed in Shanghai, but failed to appear due to Pearl Harbor). S. S. 
Hvanc, The turbulent jet, Master’s Thesis (1943), National Tsing Hua University. 

[4] N. Hu, Velocity and temperature distributions in turbulent wakes behind an infinite cylinder and 
behind a body of revolution, Chin. Jour. of Phys., (1941) (printed in Shanghai, but failed to appear). 

[5] N. Hu, The turbulent flow along a semi-infinite plate (unpublished). 

[6] Y. C. Hsten, Turbulent couette motion and the turbulent flow between two concentric rotating 
cylinders, Master’s Thesis (1942), National Tsing Hua University. 

[7] P. Y. Cuou, Theory of isotropic correlations and the turbulent flow in a windstream, Sci. Rep. of 
National Tsing Hua University (1941) (printed in Shanghai but failed to appear). 

[8] Tu. von KARMAN and L. Howartn, Proc. Roy. Soc. Lond. (A), 164 (1938), 192-215. 

[9] P. Y. Cuou, Pressure flow of a turbulent fluid between two parallel infinite planes (to appear 
subsequently in the Quarterly of Applied Mathematics). 

[10] G. I. Taytor, Proc. Roy. Soc. Lond. (A), 151 (1935), 421-478. 

[11] G. I. Taytor, Jour. of Aero. Sci., 4 (1937), 311-315. 

[12] H. L. Drypen, G. B. ScuuBaver, W. C. Mock, Jr. and H. K. Skramstap, N.A.C.A. Report, 
No. 581 (1937). 

[13] Pranptt and ReicHarpt, Deutsche Forschung (1934), 110; analysis by G. I. Taylor, Proc. 
Roy. Soc. Lond. (A), 151 (1935), 456. 

[14] G.I. Taytor, Proc. Roy. Soc. Lond. (A), 157 (1936), 537-546. 


55 


QUANTITATIVE INTERPRETATION OF MAPS OF MAGNETIC 
AND GRAVITATIONAL ANOMALIES BY 
MATHEMATICAL METHODS* 


BY 


E. G. KOGBETLIANTZ 
Lehigh University** 


1. Introduction. In geophysical prospecting for oil and other minerals gravita- 
tional and magnetic anomalies corresponding to geological phenomena are mapped. 
The problem of the quantitative interpretation of such empirical maps consists in 
the determination of numerical values for all the geological parameters (depth, thick- 
ness, slope, density, intensity and direction of magnetization, etc.) which characterize 
a tectonic structure or an ore-body. To illustrate the possibility of such an interpreta- 
tion it is preferable to avoid the complications involved in the mathematical study of 
maps of complex anomalies. The complex anomalies are due to the coexistence in 
the same region of many different geological phenomena. Resulting from the super- 
position of many simple anomalies, they can be resolved into their simple compo- 
nents each of which corresponds to a single ore-body or tectonic structure. This 
resolution is the first step which must be performed, since no interpretation of a com- 
plex anomaly map as such is possible. The problem of resolution is a very important 
one since in most cases we have to deal with complex anomalies, the simple ones being 
exceptions. Special methods devised by the author solve this important problem, but 
they are not discussed in this paper which deals with the quantitative interpretation 
of a simple anomaly map. We study here two cases: an axial anomaly created by an 
anticline and a centered anomaly corresponding to a salt dome. They are sufficiently 
simple and at the same time have great practical importance. 

A new method of interpretation, based as the usual methods on the theory of po- 
tential but essentially different from them, is introduced in this paper. In the usual 
methods! systematic use is made of individual values such as maxima, minima, 
zeros, inflection points, etc., of the observed and plotted quantity as well as of their 
distances. The use of such remarkable values and distances is founded on the tacit 
assumption that they reflect exclusively the physical action of the unknown structure 
or ore-body whose study is the object of the interpretation. This assumption is per- 
missible for the anomalies of large magnitude but it is doubtful for those of average 
magnitude and completely wrong for small anomalies. 

In the past, geophysical prospecting by gravitational and magnetic methods was 
directed mostly toward the study of important, clearly pronounced anomalies of 
large magnitude which correspond to more shallow deposits or to big well defined 
tectonic structures. But now the geophysicists are obliged to deal with more difficult 


* Received April 11, 1944. 

** On leave of absence. Now at The New School, New York. 

! For examples of these usual methods see L. L. Nettleton, Geophysical prospecting for oil, McGraw- 
Hill Book Co., New York, 1940, Chapters 6 and 12, and also H. Shaw, Interpretation of gravitational 
anomalies, Trans. Amer. Inst. Min. Met. Eng. 97, 271-366 (1932). 
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cases and they have to interpret maps of small anomalies. Thus, we must study the 
obstacles which naturally lead to an erroneous interpretation of an anomaly map by 
the usual methods of remarkable values and distances if the magnitude of the anomaly 
is not very large. 

2. Punctual anomaly. An anomaly map is plotted on the basis of measurements 
made at isolated stations. It is supposed to be generally correct, and the usual correc- 
tions required by the topography of the region, the fact that the earth is not a sphere, 
etc., are supposed to have already been made. We must emphasize that a quantita- 
tive interpretation presupposes an accurate correction for the so-called regional anom- 
aly, and it cannot be expected to yield good results if the latter correction is made, 
as is customarily (and very unfortunately) done, simply by smoothing arbitrarily the 
experimental curves. Special methods exist which ensure a very accurate correction 
for regional anomaly by deducing it from the map itself, but this important question 
cannot be discussed here. All usual corrections having been made, each individual 
value obtained at a station is the combined effect of two anomalies: 1) the anomaly 
caused by the tectonic structure or the ore-body the study of which is the purpose 
of the interpretation, and 2) the anomaly generated by local irregularities of mass 
distribution or of magnetization intensity in the immediate vicinity of and under the 
point of measurement. This strictly local anomaly—we propose to call it “punctual 
anomaly”—is in general very small. It affects only a small area around the point and 
it is precisely this punctual anomaly which is responsible for perceptible variations 
in the value of the observed quantity which occur for small displacements of the ap- 
paratus used around the station. In fact the apparatus used now are extremely sensi- 
tive and we cannot neglect any more the existence of punctual anomalies. There is no 
correction at all for them since the punctual anomalies, affecting every observed in- 
dividual value, cannot be evaluated. If the magnitude of the studied anomaly is 
large, the punctual anomalies are negligible and the map can be interpreted with the 
aid of remarkable values and distances. The positive results achieved by the old 
interpretation methods must be explained in this way. But, if the magnitude of the 
anomaly is small, the punctual anomalies not only modify the extremal values but 
they also displace them, altering all the distances used in the usual interpretation 
methods. Since these old methods express all geological parameters in terms of re- 
markable values, their distances and ratios, it is plain that punctual anomalies render 
these methods completely useless in the interpretation of small anomalies. It is a very 
important though often disregarded fact that the interpretation based on isolated 
values can in general be only qualitative and gives exactly nothing in case of small 
anomalies. This important fact explains the lack of success in dealing with maps of 
small anomalies and is the reason for the actual ineffectiveness of geophysical pros- 
pecting in discovering new oilfields in U.S.A. New methods of interpretation well 
adapted to small anomalies are now necessary. They must be introduced into practice 
if the geophysical prospecting by gravitational and magnetic methods is to be applied 
in the future. 

The interpretation errors caused by punctual anomalies can and must be elimi- 
nated and there is only one possible way to do it. Considered together, the punctual 
anomalies in the region covered by the measurements have a random distribution; 
they oscillate about zero and are independent one from another. Consequently, they 
must undergo an almost total compensation if we form the average value of some 
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function of the mapped quantity, the average value with respect to the whole map, 
using all the observed values at a time. We can eliminate the harmful influence of 
punctual anomalies and compensate at the same time for the possible residual ob- 
servation errors only by combining all observed values in an integral. For any map 
in general the average values express much better the action of the phenomenon under 
study than do the individual values and their distances. In other words, interpretation 
rules and formulae based on the systematic use of integrals give much more correct 
quantitative results in all cases and for all maps. Rules based on special individual 
values hold only in very rare and exceptional cases of big structures such as, for in- 
stance, the shallow salt domes of Texas, the Kursk and Kirunawaara iron ore-bodies 
or the Great Rhodesian Dyke. 

The method described in this paper uses exclusively the average values and, in 
particular, moment functions and moments of the observed quantity and of its square. 
This method was applied by the author in France and in Iran with good results. The 
cases studied here are chosen only for the sake of brevity. The method is elaborated 
for the most general cases of complex anomaly maps obtained as result of magnetic 
or gravimetric survey of a completely unexplored region. 

3. Center of gravity and first moments. The problem of locating the center of 
gravity C of disturbing masses is a fundamental one, and its solution is the first step 
of every interpretation. We solve it in the general cases of an axial anomaly and of a 
centered anomaly. The coordinates x*, y*, 2* of C are expressed with the aid of mo- 
ments of the observed quantity Q, that is in terms of 


Om =f eoaz 


for an axial anomaly and in terms of 


Onn = ff x™y"O(x, y)dS, Om = ff r™O(r, o)dS, (r? = x? + y?, x tang = y) 
P Pp 


for a centered anomaly, the double integration being extended over the infinite plane 
xOy denoted by P. 

Axial anomaly. If the geologic feature being considered is much longer in one di- 
mension (strike, axis of anomaly), the corrected map of the axial anomaly created by 
such a structure is a family of nearly parallel lines and the interpretation deals with 
a curve describing the behaviour of the plotted quantity on a typical profile perpen- 
dicular to the anomaly axis. Cartesian coordinates x, y, z are introduced with the 
z-axis directed vertically downward, the origin O being at the surface and the x-axis 
being perpendicular (the y-axis being parallel) to the anomaly axis, as shown in Fig. 1. 
The excess of the density of disturbing masses over the density of their environment 
is called the density-contrast and is denoted by o (o $0). We represent the disturbing 
structure as a homogeneous cylindrical body of normal cross section S, denoting the 
area of S by A. Ata point (x, z) in the plane y =0 the potential U of the body is given by 


U(x,z) = —- eff log [(x — £)? + (2 — £)*]dS + const., (1) 


where &, ¢ are running coordinates on S and k = 2fe, f being the gravitational constant 
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(66.7 X10-® cgs.) and o the density-contrast. On the surface (plane P), z=0, the par- 
tial derivatives U, and U,=Dg of U are given by 


U,— iDg = — ef f (x — p)— dS = bg log (x — p)df, (2) 


where p=£+7¢ and I is the boundary of S. The second derivatives U,, and U;, are 
called the curvature K and the gradient G respectively: K = Uzz, G=U;,. Thus 


K-—-iaG= ef f (x — p)~*dS = eg (x — p)— df. (3) 
s r 


Using in (2), (3) the binomial expansion, we deduce for large | x| the approximations 
which hold for any form of the section S. Their first terms for instance are 


U,~ — kAx, Dg ~ kAz*x~, G ~ — 2kAz*x-3, K ~ kAx, (4) 
where 2* is the depth of the center of gravity C. Denoting an arbitrarily chosen origin 


of the coordinate x’ on the profile by O, we regard the function Dg = Dg(x’) as known 
from the measurements. To find x* =O0O*, z2* =O*C we shall use the first three mo- 


ments of Dg, . 
Mo -f Dg(x’)dx' = rkA, M, -{ x'Dg(x')dx’ = rkAx*, (5) 


J tee) — kAz*|dx = rh ff (g? — ¢?)dS, (6) 


where in (6) the origin of the coordinate x is the point O*, the projection of C on the 
profile, x* being considered as already found with the aid of (5). In fact, from (2) we 


deduce that 
x’U,+ kA — ix’Dg(x’) = ef f — x’)"pdS. (7) 
s 


Integrating (2) and (7) with respect to x’ in (— ©, ©) and observing that the in- 
tegral of (x’—p)—'dx’ equals im since [ in p=§+7£ is positive, we have (5). To prove 
(6) we integrate in (— ©, ~) for x*=0 the relation 


a(xU, + RA) — i(x*Dg(x) — kAz*) = ef f (p — x)—'p’dS 
s 


and compare the coefficients of imaginary terms. 

From (5) we deduce the rule x* = _M,/Mpo. However, in practice the integration 
can be carried out only over a finite interval (—R, R), where the known length R is 
at least four or five times the depth 2* of C. The contributions from the intervals 
(— oo, —R), (R, ©) can be computed by means of the expansion of Dg(x) in powers 


of x, 


Dg(x) = kAs*x-*{1 + 2ci2*a-! + (321 — Cos)2*2x-* + O(x-*) j, (8) 


where the constants C. are given by 


feymtean = ff enras 
8 
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In general, unless S is very irregular, we have 2*cn+x*, 2**ca.=x*?, cos= 1, whence (8) 
takes the form 
Dg(x) = kAs*x-?{1 + 2x*a! + (3a*? — 2%)? 4 O(x-*)}. (9) 


If M#* and M* denote moments computed for the interval (—R, R), i.e., 


R R 
Mt = Dg(x’)dx’, Mt -f x’Dg(x')dx’, 





—R R 
and if we neglect terms of the relative order O{(s*/R)?} , we find that 
22* 42* 
wkA = My = (14+—) mp, M,= 1+—) M*, (10) 
aR mR 
whence 
» (1 + =) Ms (11) 
z* = — . 
aR Mt 


M¢ and M* can be obtained from the experimental curve for Dg(x’) by means of a 
planimeter. From (11), we note that M*/M¢ is a first approximation for x*. 

If we set x*=0, by (9) we easily see that the contribution from the intervals 
(—o, —R), (R, ©) to the integral on the left side of (6) is of order (2*/R)?. If we 
neglect terms of order (2*/R)?, (6) can be written in the form 


—Me=M += ff (@ — ¢)aS 
tested 9 FM 


where, in the usual notation, 


R 
Ms -{ x*Dg(x)dx. 


—R 


If we consider only those cases in which the horizontal dimension of S is much smaller 
than its depth, then ff,(g?—{*)dS+ —2**A. Using this, and substituting for Mo from 


(10), we have 


M+} (x? — 4)z* 
ail! neh | (12) 9 9" 
2RM¢ 2rR 








~ _ 











Equations (11) and (12) permit (x,2) 
us to compute x* and 2* by succes- 2 
sive approximations, the successive 
values being denoted by x,*, 2," fi 
(n=1, 2,3, - - - ) and the correspond- «e,6) 
ing positions of O* (Fig. 1) by O,*. is | 
The steps are as follows: (a) We 
choose a position for O and compute: 
M,*, M3; integrating over the inter- 
val —R>x’<R with a planimeter. 
(b) We obtain x* by setting 2* =0 in 
(11) and plot O¥. (c) We compute M* Fic. 1. 











N 
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and new values of M#* and M*. (d) Using (12) with z*=0 in the right member we 
compute 2;*. (e) Using the values of M,*, Mi*, Ms obtained in (c), and replacing 2* 
in the right members of (11) and (12) by 2*, we obtain x*, 2* and plot O.*. The steps 
(c) and (e) are repeated until a stabilization of the points 0,* and the values x,*, 2,* is 
reached. 

The same method can be applied to maps obtained by means of a torsion-balance, 
giving the curves of the gradient G(x) and the curvature K(x). Denoting the mth 
moment of G by H,, we have Hy=0. Also, integration by parts yields H,;= — Mo, 
H,= —2M;, H3=—3Mb2. For the corresponding reduced moments we have Mo = 
—A*(1+22*27—R-!), 2M* = —H#(1+22*/rR), 3M = —HX+2kRAR2*. Thus (11), 


(12) can be transformed into 


: Ht 22* r hH# ps* 
o's = 1+ — |}, , oa (13) 
2H; wR RH; R 


where A =32/(39r—1) =0.187, uw =3(30? —242r+8)/(30?—7) =0.485. Equations (13) 
permit us to determine the center of gravity C from the gradient map only. 

The first two moments Lo and L; of the curvature K vanish. Since for | x| very 
large K(x)~kAx~*, we define the moment Lz and the reduced moment L* by the 


oo Hy, R H* 
LI -f x*K — —  }dz, L¥ -f “kK — —) dx. 
= wT / = TT 


If we multiply (3) by x?, subtract kA from both sides and integrate over (— ©, ~), 
we find that Lz=2//:2*, the contribution from the intervals (— «©, —R), (R, ©) being 
—6kAz*?/R. Since H, = H*(1+42*/7R), we then have for 2* 


H#z* = aL#(1 + Bz/*R), (14) 





integrals 


where a =41?/(7?—4) =0.84, B=2/(r?—4) =0.535. Equation (14) supplies a control 
on Eqs. (13). 

The magnetic anomaly created by a cylindrical body of section S is related to the 
gravitational anomaly generated by the same body, and the equations relating to 
maps of G and K can be transformed into equations relating to maps of the horizontal 
and vertical components X and Z of the abnormal magnetic field created by the 
body. If J and w denote the magnitude and inclination of the magnetization vector, 
we have the classical relation (Poisson) 

k(X + iZ) = 21(K + iG)e~**, (15) 
where k =2fe, f being the constant of gravitation and o the density-contrast. Multi- 
plying (15) by x and integrating over the interval (— R, R), we obtain for the reduced 
moments X#*, Z#* the relation k(X*+7Z*) =21(L*+iH*)e—*. Since Li =0, we easily 
find that rRL#* =4x«*H}#*(1+42*/7R). Thus 


4x* i. 
kX* = 21H* (sin y+ "9 cos v), kZ* = 21H*# | cosy — —. sin v), 
1 


qT 
whence we obtain for the two parameters JA and y, 
ITH} 1 


cx* 
fw = — ow — (Xf? + Z7*)*, Z* tany = xe(1 ox =), 
kr 2r 
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where c is defined by the relation tX*#Zi*c=4(X#2+Z}?). To find x*, 2* we need 
second moments. Since the principal term of the right side of (15) involves K, and 
K~kAx~, for large |x| we have rX ~Zj*x~?, rZ~ — Xx. Therefore we define the 
second reduced moments X¥,Z* by the integrals 


R R 
X3 = f (x?2X — Zf*a-')dzx, Zi = (x°Z + Xfa)dx. 
ad —R 

Their numerical values can be obtained from the experimental data as areas under 
curves deduced from the curves X = X(x), Z=Z(x). 

Integrating (15) after multiplication by x*, and using the defintions of X#*, Z#, Li, 
we obtain 

k(X# + iZ#) = 2e-¥ [LF + i(HF — 2RL Fr )]. 

Substituting in this result the values of H#, L# obtained by solving (13) and (14), 
and using the relation 2rRL* = 8x*H#*(1+42*/rR), we find that 


s* =a(l+ys*R)N.,  2* = a(1 + Bs*RN, + 4x*4 RA, 
where @ and £ are as in (12), y=2(16—2?)/(2? —47) = 0.665, and 
Nz = (XPX4 + Z7Z5)(X + Zi), N, = (XSZ* — X2Z*)(X"? + 27)". 


The numbers N,, N, can be deduced from the maps. Thus the moments of X, Z give 
the four parameters JA, y, x*, 2*. } 

It is to be noted that the above results can be obtained with the aid of the theory 
of Fourier transforms. All our results are based on (2), which can be written as a 
Fourier transform. Now p=£+2¢, and since min.(¢) >0, we have 


rc ) 0 
i(op - x= f eta) dt = f e~ #2) df, 
0 0 


This proves that Dg+iU, is the Fourier transform of a function w(t), vanishing for 
positive ¢ and defined in the interval (— ©, 0) by the relation w(t) =k(2r)"?f,e-# ‘dS. 
On the other hand, Dg—iU, is the transform of w(—#), which vanishes for negative ?, 
and Dg(x) appears as the transform of the function f(t) defined for all values of ¢ by 





f(t) = an ff efft-f tlds; (16) 
s 
Dg(x) = (on) f e***f(2)dt. (17) 


This expression will enable us to find easily the moments of Dg. The moments of the 
square of Dg, which will be required presently can also be deduced easily from (17) 
with the aid of the Parseval theorem. 

Centered anomalies. In the case of a centered anomaly, we represent the disturb- 
ing structure as a homogeneous irregular body B. Cartesian coordinates x, y, 2 are 
introduced, with the z-axis directed vertically downward, the origin O being arbi- 
trarily chosen on the plane P (Fig. 2). C(x*, y*, 2*) is the center of mass of B, and O* 
is its projection on the plane P of the map. 

The gravitational anomaly is Dg(x, y). In the old methods, O* is placed at the 
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maximum of Dg. If the body is a solid of revolution with a vertical axis, this is cor- 
rect; but if the body is irregular or inclined, the maximum of Dg occurs somewhere 
above its uppermost part. In the present paper, we shall locate C by means of in- 
tegrals. 








dV 
(87,5) 








If 5V is an element of volume of B at a point (£, 7, ¢), and if 5[Dg(x, y)] is the 
contribution to Dg from 6V, then 5[Dg]=fot6V[(x—£?+(y—n)?+@]>. Now 
J fri[Dg|dS =2xfosV, which is independent of £, n, ¢. Integration of this over B 
yields [{pDgdS=arkV, where k=2fo and V is the volume of B. This result holds for 
any homogeneous irregular body or bodies. 


Because of symmetry //, [x—E+4(y—n) ]6[Dg(x, y) ]}dS=0. Thus 


ff (x + iy)6[Dg|dS = (& + in) ff sslas = wk(E + in)iV, 
e P 


and integration over the body B yields 


ff DgdS -ff xDegdS, x ff Deas = f f yDegdS. (18) 
F P Pp 


These are two equations for x* and y*. They hold for complex anomalies as well as 
simple ones. In practice, integration can be carried out only over a finite part of the 
plane P. We choose that part lying inside a circle with center O and radius R, where R 
is a constant at least four or five times the depth 2* of C. The equations of this circle 
are r=R, z=0, where r?=x?+y?. We denote its interior by L. The contribution to 
the above integrals from the infinite region r= R can be easily evaluated, since at such 
large distances the gravitational action of the body is approximately the same as that 
of a punctual mass a V located at the point C(x*, y*, 2*). Therefore, for r2R we use 
the approximate formula 


Dg = foz*Vr-*{1 + 3(x* cos @ + y* sin 6)r—! + (9x*? + 9y** — 62**)(2r)-? 
+ 15[2x*y* sin 20 + (x*? — y**) cos 20 |(2r)-? + O(r-*) . (19) 


where 7, 6 are polar coordinates in the plane P, with origin at O. Neglecting terms of 
order (z*/R)- and higher, we have with the aid of (19), 
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(1 = =) f J Deas = f J reas, 
(1 = Vf fie + iy)DgdS = fc + iy)DgdS. 


Thus (18) can be written in the form 


s* 3g** 
(x* by* DgedS = {1 — ff iy) DedS. 20 
+ ix) ff peas = (14+ 5+) ff cet inne (20) 


Equations (20), applied in the first approximation with 2* =0, give a first position O; 
for O*. If we choose the origin O at this point, then d=O,0* is small. 

To obtain an equation for z*, we integrate r?DgdS over L. Neglecting terms of 
relative order (z*/R)*? and higher, we need only the two first principal terms of this 
integral. Thus we can evaluate the contribution of an elemental volume 6V at (&, 9, £) 
by integrating r?6[Dg]dS over the region r’ SR instead of L, the origin of polar co- 
ordinates (r’, 6’) being at the point (£, 7) above the point (&, 7, ¢). In fact, the differ- 
ence between two integrals over Z and r’SR is of relative order (2*/R)?. Now 
r?=r'2-+-9?—2pr’ cos (@—80’), and integration over r’<R gives 





32* 
2R 





f : r°6[Dg|dS = rkR6V {¢ + (p? — 2¢2)R- + £0(e*2/R)}. 
L 


Integrating this result over V and replacing the integral of the second term by its 
approximate value —2kmr2*?V, we obtain 


f f r°DedS = rkVz*R{1 — 22*R-! + 0(2**/R)}. 
L 


Dividing by R/{:DgdS =xk VR(1—2*R-'), we find that 


or ff Deas = (1 + 2*R-") Jf roses . (21) 


The term 32*2R-? must be added to the factor 1+2*R- if the terms neglected are of 
order (2*/R)* and higher. 

When the measurements are performed with the aid of a gradiometer or torsion 
balance, the resulting maps of U,, and U,, give not only x*, y*, 2* but also a control, 
since each of these two maps can be used to locate the point C. Applying the same 
reasoning as for Dg, i.e., first integrating 6U,, and 6U,, corresponding to an elemental 
volume 6V, and then integrating the result over V, we find that 


1 
off xU dS = —ff x?U dS, * ff yU,dS ff xyU,yd5S, 
L 2 L L L 
1 
ff xU dS ff xyU ,dS, ff yUydS —ff yy dS. 
L L L 2 L 


All the integrals in (22) have the same reduction factor 1+32*/2R. Hence it does not 
appear. The third reduced moments give equations for 2*. They are 


(22) 
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73* 72* 
ore ff sU 4S = 8(1 + Sf x8U 2dS 24(1 + ‘ \ff xy7U dS, 
L R L R L a 
a (23) 
one ff yUydS 8(1 + Sf y?Uy.dS a(1 + 
L 6R - 





ll 











72* 
x*yU, aS. 
SE i 


If we neglect only terms of order (2*/R)* and higher, the term 112*?/(18R?) must be 
added to the reduction factor 1+72*(6R)—'. Equations (22) and (23) are general. For 
a solid of revolution with a vertical axis, they reduce to 


: 3c2* ; 
(x* + iy ff GdS -( on Sf (x + iy)GdS, (22*) 
L 
2* 112*? 
3Rot f rGdS = 2(1 +4 sa LS. riGdS, (23*) 
L 


where ¢ is such that Ref {,GdS=Jf,rGdS. 
In the general case of a complex magnetic anomaly we assume that the magnetiza- 


tion vector is the same for all particular isolated bodies which create this anomaly. 
Its intensity J, inclination y, and azimuth ¢ are three unknowns; J cannot be sepa- 
rated from the total volume V and it is the product VJ which is deduced from the 
maps; ¢ is defined with respect to arbitrary cartesian axes of x and y on the surface 
of the earth; X, Y and Z are the components of the anomaly. 

We shall now deduce expressions for the six parameters VI, y, ¢, x*, y*, 2* which 
characterize the magnetization vector and locate the common center of gravity of all 
the disturbing magnetic bodies. This can be done by computing the moments of 
X, Y, Z from maps showing the distribution of X, Y, Z over the surface of the earth, 
and by use of the classical Eétvés formulae, 


fo(iX + jY + kZ) 








re) 0 0 
= 1(i —+j—+ «—) [((U.cos¢ + U,sin¢) cosy + U,siny], (15*) 
Ox fe) Oz 


where i, j7, & are unit vectors on the coordinate axes. We shall use moments of the 
first four orders, denoting them by subscripts. For example, Xmn=JffrXx™y"dS. It 
is to be noted that X,,, are reduced moments. Neglecting terms of relative order 
(z*/R)?, and using the same method as in the case of (21) in the computation of in- 
tegrals of the second derivatives of U, we obtain 


— Xm sec ¥ sec d = — Yoo sec y csc d = 2Zy9 csc YW = 2nVI/R. 
The first approximate values of VI, ¥, @ are then given by 
-1/2 1/ 
tan 6 ="Feo/Xse, tay = BZe(Xec+ Fu)» sVI = R(X + YootiZe) - (24) 


The zero moments are small quantities, and approach zero as R approaches in- 
finity. Hence (24) yield poor approximations. Nevertheless, the first value of ¢ must 
be used to rotate the axes of coordinates, directing the x-axis nearly parallel to the 
horizontal component of the magnetization vector, so that ¢ will be very small. 

Among the six moments of the first order, we do not use Xo and Yio since they are 
of order 1/R compared with the other four, which are given by 
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Fass - ~ sin y(1-+ = cot cos#), ak. = — sin y(1-+ cot sin 6), 
2nVI 2R 2nVI 2R 
410 = — cosy (cos - se tan v). on =— cos (sin o- x tan v). 
2xVI R 2nVI R 
Hence we obtain better values for VJ and tan y, 
QeVI = (Zin + Zu + X10 + ivan(1 + ~), 
tany = ho? + You) (Zi +2) "(1 te =), 


where 
8c, = 3 sin 2p(x* cos ¢ + y* sin d) = 6 cos? y-cs 
4co = (cot ¥ + 4 tan y)(x* cos ¢ + y* sin d) = (4 + cotan? )cz, 


ll 


—3/2,_ 2 2.1/2.2 2.~2 2 2 2 
Cés=2 (Xiot Yio) Zor +210) [(Z:20 — Zo2)(Zo1 — Z10) — 4ZuZu:Z10). 
The coordinates x*, y* of the point O* are obtained with the aid of second mo- 
ments, as indicated by the equations 


%* 2 


2 = c y 
3(Z10 + Zn) ‘[(Ze0 — Zo2)Z10 + 2ZnZ1| (1 te “) + R tan y cos ¢, 


* 


x 


(25) 


7 


all 2.-1, c3 ies P 
$(Z10 + Zo1) [2Z02Z21 — (Zoo — Z02)Z01] i- R + PR tan y sin ¢. 


The expression 





. 5 , ‘ (1 =) 
oS == cos sin _ 
aa ial om? 5R 


shows that Yoo is zero if the x-axis is parallel to the horizontal component of the 
magnetization vector J. If Y2o is different from zero, the value of ¢ is given with good 
precision by the important equation X oz = Yeo tan ¢. Thus, by use of this equation and 
Eqs. (25), we can change our axes so that x*=y*=¢=0. Computing the third mo- 
ments in this new system of coordinates, we find that Xn = Xe = Yso= Yu =Zn =Zos 
=0, 3Z12=Z30 = X30 cot yy, 


3. = , = xX = ‘ — << = 2* a a! sin Ww 
12 21 30 03 R 


These expressions give for tan y the eight values X30/Zs0, Yn/Zu, etc., and the mean 
of these eight values can be considered as the final value of tan y. Now, with this 
choice of axes, the second moments take the values Xi = Yoo= Yoa=Zu=0, 


: 5 162* ¥ 5 R(1 —) y 
d to = R 1 —> ’ 2 = oe = cos ’ 
7 ( SR ) ts Td SR 


3 82* 142*\ | 
Yu = R(1 - =) cos Y, Z20 = Zoe ia sr(1 = SR sin y. 
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Dividing any third moment by one of the five second moments different from zero, 
we obtain an approximate value for 2*, and the mean of all such values can be consid- 
ered as the final value of z*. Thus, we have not only solved the problem, but have 
also found a control of the solution, since the degree of concordance of many different 
values found for the same parameter, such as 2*, characterizes the reliability of the 


solution. 
Fom a practical point of view, a solution involving only the moments of the verti- 


cal component Z is important, since in many cases magnetic surveys are limited to 
measurements of Z only. In such cases we can find Wy and @¢ and locate the point O* 
by use of Eqs. (25) together with 


y* — x* 
Zy) tan gd = Zu( 1 oa nae tan v), 


2.1/2 C 2 241/2 
SRZi0 + Zo) tany = (1 + Ie — 2x*Zi) + Zor — 2y*Zur) | 


where ¢c, is such that 5Re,= 142* —(x* cos ¢+* sin d) tan y. For VI we have 


1/2 


C4 2 2 
107RVI = (1 + «ein + Zo2) 


Also, 2* is given by 





xn 
| 


8 72*\ Z30 2 2.1/2 
*= ->(1+ \F Gi + 2h) ’ 
6R/R 


where the moments are related to the origin O* with ¢=0. It is plain that the accuracy 
of the approximate values given by the above equations increases with R. 

4. Interpretation of the map of an anticline with the aid of moments. Let us repre- 
sent the structure of an anticline as a cylindrical body of normal cross section S. We 
have first to choose for S some simple geometrical form which expresses the gravita- 
tional action of the anticline adequately. Our choice is the region between two con- 
centric elliptic cylinders, the major axes of which have an angle a of inclination 
(Fig. 3). This choice is good, except in the vicinity of the deepest part of the region, 
which part is at a considerable depth; the author has verified that such a choice leads 
to good results in practise. We shall denote the outer and inner ellipses by £; and Ey, 
respectively, and their semi-axes by a, b and ya, yb, respectively, where ¥ is a positive 
constant less than unity. 

Equations (10)—(14) give values for kA, x*, 2*. We shall now show how y, @ and 
the eccentricity e can be deduced from gravitaticnal anomaly maps. 

Let us replace e by a parameter r = 22*/c, where 2c is the focal distance of the outer 
ellipse. We define a function n(s) by the relation 


n(s) = My M(—s)D, (26) 
where s is a constant and 0<s<1, Mp is as defined in §3, M(—s) is the moment of Dg 


of order —s, and D is the zero moment of [Dg(x) ]?: 


M(— s) -f | «|-*Dedx, D -f [Dg(x) ]*dx. 
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The graphs of |x|~*Dg(x) and [Dg(x) ]* can be plotted from the experimental curve 
of Dg(x), and thus the reduced moments can be computed with the aid of a planime- 
ter. On the other hand n(s) can be tabulated for various values of s. We shall require 
n(4), n(4), n(2). 

Applying the theory of Fourier transforms to (16), (17), we have 


f Dg(x)e*"dx = rk ff etlel+ittgs, (27) 
S 


—2 . 


2f [Dg(x) }2e-‘"dx = re ff ff asas' f ett (du, (28) 
ry Ss Ss 2 

















Oo O a 
SAAN CAAA DAK N PE SNUNENSNS SUNT NN OQ 
B 
p=2° A 
a 

A B’ 

x" 
& 
Fic. 3. 


where $(u) = &+(£’—&)u —{'| u| —{| t—u| , and (&, ¢), (&’, ¢’) are two sets of running 
coordinates on S. We note that 


f e*u(y — t)-8dy = (iv) e'T(1 — 8) 
t 


for £20, 0<8<1, Reliv] <0, where Re[iv| denotes the real part of iv. Consequently, 
the application to (27) of fractional integration of order 1—64 yields for ¢20, 


f Dg(x)| x|*e#dx = rk f f p>e”'dS. (29) 
— Ss 


The inversion of order of integration is permissible, since all integrals are absolutely 
convergent. Replacing in (29) the exponent 6—1 by —s, (0<s<1), multiplying both 
sides of this by e***/? and letting ¢ tend to zero, we obtain 


M(-— s) -{ | x 


—o 


—*Dg(x)dx = rk sec (}4s)Re[(A; — s)], 





where (A; —s) =eirsl2 ff p-dS. For our region S we have 
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®(A; — s) = A(s*)-*[F(4s, 3 + 45, 2; — 4r-%e-*a) 
— y°F(4s, § + 4s, 2; — 4r-*y%e-***)], 
where F(a, 6, c; x) is the classical hypergeometric function. The function M(-—s), 


which depends on the four parameters s, a, y, 7, can be tabulated for —1<s<1, 
0s7 31, 0SaSjr, r22 sin a. We shall require tabulations for s =}, 3, 2. 


If in the right side of (28) we carry out the integration with respect to u for 20, 
and then let ¢ approach zero, we find that 


2D = 2f [Dg(x) ]*¢x = ik? ff ff (p — p’)—dSdS’, 


where p=£&+2¢, p’ =&’—i¢’. In the present case, we have 
D = Mi(s*) ¢(a, 757), 
where 
(a, ¥; 7) = fle, e*; 7) — Rel flye-, e; 2) ] + flye-*, ve; 1), 
3f(u, 0; 7) = 16ru*v®(u + v)-?[C\E — C2K + C3KE(b, ’)], 
K and E being complete elliptic integrals of the first and second kind of modulus A 


given by 4uv =)?[r?+ (u+v)?], X’ and b being such that A’ = (1—A?*)"/2, 2(uv)"/2dn(b, Xd’) 
=i-(u+v), and Ci, C2, C3 such that 

C3 3n’[1 — d’sn(b, d’)cn(b, d’)dn(b, d’) a 

C, = 3dn?(b, \’) + 3d’sn(b, ’) — 1 — 0X? + OC, 

C2 = 3X/2dn?(b, d’)cn?(b, N’) + 3d/sn(b, dX’) [1 + d/2cn?(d, dN’) ] + 2 + OC. 


Thus (a, y; 7) can be tabulated. 
We have now derived expressions for M(—s) and D occurring in the right side 


of (26). Substitution yields the equation 
E(s; a, y; 7) = n(s) (30) 


where E(s; a, y; 7) =[(a, ¥; 7) ]-*m(—s; @, y; 7), with 
m(— s;a, y; 7) = sec (44s)A—!z**Re[®(A; — s)]. 


In (30) we set s=}, 3, 3, to obtain three equations which we can solve for the re- 
quired quantities a, ¥, r. 
Once a, y, 7 have been determined, we can deduce a new value for 2* from the 
relation 
‘ 2._—l1 
z*= MD (a, 7; 7). 


This serves as a control on the value obtained from (12). 

In actual computations based on experimental data, we can use only a finite re- 
gion on the surface of the earth. Hence reduced moments must be introduced, as 
before, and the parameters of the problem must be deduced by successive approxima- 


tions. 
The same procedure can be used in treating maps of the gradient G and the curva- 
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ture K. Each map leads to an independent set of values for x*, 2*, a, y, 7, RA. In the 
case of a magnetic anomaly, the map of X or of Z could be used. In all cases the mo- 
ments can be expressed in terms of ®(A; —s). 

We note that the above method yields kA. Now A is the area of the cross section 
and k=2e, where f is the gravitational constant and o is the density-contrast. We 
define the average thickness T of the disturbing layer as the geometric mean of the 
extreme thicknesses; JT =(1—v)(ab)!/2. Since A =2ab(1—y?), we then have A(1i—vy) 
=m(1+~7)T?. Hence, if ¢ is known, A can be found and then T. Often T is known. 
In this case, A can be found, and then a. 

5. Interpretation of a centered anomaly created by a salt dome. Before applying 
the new method it is interesting to see what can be obtained by the old method. 
Some results can be achieved if we disre- rey 
gard the cap-rock, neglect the slope of the 
flanks and omit the depth of the salt dome 
base considering this structure as a homo- Pp 
geneous vertical infinite circular cylinder 
whose top is at the depth p (Fig. 4). The 
interpretation problem is reduced to find- ae 
ing the three parameters: k(o=k/2f), p 
and the radius a of the dome. Instead of a 
we consider the angle w, letting a=p tan w. 
With the origin at O*, just above the cen- 
ter of the top circle (maximum of Dg, zero 
for the gradient G) the expressions for G 
and Dg are: 

G(r) = k(a/r)'/?[2X-1(K —E) — \K], 
Dg(r) = k[(C:E + C.K — xpf(r)], 


where C, = bp sign(a — r) + 2(ar)*d-), 

C,=\(a? —r?) (4ar)-/2?+-p[E(b, d’) —b]sign(a—r) and 2f(r) =1+sign(a—r), sign 0 be- 
ing defined by sign 0=0. The elliptic function E(b, \’) and the complete integrals K, 
E have the moduli \’=(1—A?)"? and X is defined by \[p?+(a+r)?]"/? =2(ar)"/2, 
r being the distance from the origin. The argument b (0SdSK’) is that in 
dn(b, d’)-[p2+(a+r)?]"2=a+r. Either of two curves G=G(r), Dg=Dg(r) can be 
deduced from the other by graphical differentiation or integration. Therefore, we can 
use both of them for the interpretation. 

Here it is the maximum of the product F(r) =r'l2| G(r)| which it is important to 
locate and therefore the curve G must be transformed into the graph of the function 
F(r). In fact, the equation F’(r)=0 reduces to \*E(p?+a*—r?)=0 and the maxi- 
mum of F(r) corresponds to r=ro=(p?+a*)"/*. The corresponding value of the 
modulus \ is Ao=[1—tan® (4/4—w/2)]"/2. The value of the maximum itself is 
max F(r)=Fo=k(p tan w)"/?[2\¢'(Eo—Ko)+AoKo]. Thus, we have at our disposal 
three experimental data ro, Fy and the maximum Dgo=7kp(sec w—1) of Dg(r). First 
we find the angle w, solving the equation B(w) =n, where the function B(w) is defined 
by m sin (w/2)(1+sin w—cos w)B(w) =—2?[(1+sin w)Eo—Ko], the modulus Xo 
being a function of w only. The value of the number m is deduced from the measure- 
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ments by the rule »=1)/*Fo/Dgo, this value being simply the ratio of two maxima, 
multiplied by the square root of the observed distance ro. The function B(w) is easily 
tabulated. Knowing w, we have p=ro cosw and a=p tan w. The value of k is de- 
duced from Fy or Dgo. That is all that can be obtained, using the old method, and 
it is evident that it cannot satisfy a geologist. In practice the best value for ro is 
(S/2)*/?, where S denotes the area of the closed curve, the locus of maxima of F(r) on 
all radial profiles through the origin O*. This rather rough method of interpretation 
works if it is applied to an isolated salt dome. The advantage of the inaccurate first ap- 
proximation which it gives consists in the simplicity of computations. The interpreta- 
tion requires only the table of values of the function B(w) and it can be performed as 
fieldwork and very rapidly. 

We shall now apply the new method. It is assumed that the salt dome is a solid 
of revolution. Cartesian axes are chosen, with the z-axis directed downward along the 
axis of revolution, and the origin O* on the surface of the earth. We denote the cylin- 
drical coordinates of a general point inside the solid of revolution by (p, ¢, £), and of 
a general point outside by (r, 0, 2) (2<{min.); R is the distance between these two 


points. By means of the classical formula 
2a do —) 
—=2r e~S-2)“Jo(ru)Jo(pu)du, 
0 R 0 
where J,(t) is a Bessel function, we deduce for the potential 


U(r, 2) -{ e“*W(u)J o(ru)du, 
0 


where 


k 
W(u) = rk ff e-t"Jo(pu)pdpdt = =f et"J 1 (pu) pdf, (31) 
A u Cc 


A being the region the revolution of which generates the solid, and C being its bound- 
ary. Denoting the Hankel transform of order r by H,, so that by definition 


F(r) = H,|[f(u)] = f : f(u)J (ru)udu, 


for z=0 we have U(r, 0) =Ho[W(u)/u]. Differentiation of (31) yields similar expres- 
sions for U,=Dg(r), Uz, G, K: 


Dg(r) = Ho[W(u)], Use = Hol[uW(u)], —G = Hi[uW(u)], K = H2[uW(u)]. (32) 


The advantage presented by the expressions (32) consists in the possibility of us- 
ing the theorems of the transform theory in calculating the moments and moment 
functions used in the interpretation. Because of lack of space we can give as example 
only the general expression for the moment function of Dg which holds for any form 
of the solid of revolution, the expression applied below to the interpretation of the 
gravity map of a salt dome. There is another approach to the mathematical problem 
of computing the moments and moment functions needed in the interpretation. In- 
deed, direct integration of the explicit formulae (32) is easily performed with the aid 
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of divergent and summable integrals of Bessel functions of the general type (with a 
and 0 restricted only by a+d>-—1), 


J “eT o(ud)dt & 2 [(a + b + 1)/2]{P[ — a + 1)/2)}-4 (33) 


but we prefer to use the transform theory. A salt dome creates also a magnetic anom- 
aly and, for a solid of revolution, (15*) gives the expressions for its components on 
the ground z=0 with the aid of gravitational quantities G, K and U,,. Using the re- 
sults of Section 3, we choose the origin and the x-axis so that x* = y* =0, ¢=0. Under 
this assumption 


k(X + iY) 
kZ 


I [cos (Ke? — U,.) + 2G sin pe'*], 


; (34) 
21 (sin YU... + G cos y cos 8). 


Instead of X, Y it is convenient to use in the interpretation the radial and transversal 
components A, and A, of the horizontal anomaly (X?+ Y*)"/*. Transforming the X- 
and Y-maps into the maps of A, and A,, we can consider their moments and moment 
functions as experimental data. The corresponding expressions in terms of parame- 
ters of the problem are deduced with the aid of 


k(A, + iA,) = I [cos ¥(Ke® — U,,e-) + 2G sin y]. (35) 


With the aid of (34), (35) all the rules for the interpretation of the gravitational 
anomaly can be adapted to the interpretation of magnetic anomalies produced by a 
solid of revolution. Far from the origin O* the gravitational action of a body can be 
approximately expressed as that of a material point of the same excess-mass M=oaV 
located at the center of gravity C(0, 0, z*). Thus, for large r we have the approximate 
formulae 


De ~ fMz*r, —G~3fMetr, K~3fMr3, Un~ — fMr-, (36) 


the neglected terms being of the relative order (2*/r)*. For the magnetic quantities 
A,,Z we have, denoting the volume of the salt dome by V, (the caprock does not 
create a magnetic anomaly), 


A,~ IV. (2 cosy cos 6 — 32*r— sin y), —Z ~ IV,r-*(sin » + 32*r— cos y cos 6). (37) 


The moment functions are defined by the integrals extended over the infinite plane P 
and the formulae (36), (37) are used in computing the reduction factors introduced 
by the integration over the finite area r= R, denoted by L. 

The formula (32) Dg(r)=HoW(u) gives immediately W(u) =HoDg(r), that is, 


2nW(u) = ff Dg(r)Jo(ru)dS. (38) 


If the form of the solid is considered as known, only its dimensions being asked, i.e., 
if the expression of W(u) as function of the parameters is prescribed, the relation 
(39) becomes the source of equations since for every particular value of the arbitrary 
parameter u (which is the reciprocal of a length) it is an equation in the unknown 
parameters of the problem. In practice the numerical value of the second member is 
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obtained by integrating over L, and adding to the result the contribution of the 
infinite area r2R. This contribution, computed with the aid of (36), is equal to 
2xf M(z*/R)[Ji:(Ru)/Ru] where, as has been seen in Section 3, 2xfM=rkV is the 
value of the first moment of Dg. In this method the map of Dg is transformed, for 
each value of u, into the map of the product Dg(r)Jo(ru) and the integration over L, 
giving the second member of (38), is performed on this map. Let us consider as an ex- 
ample the special case when the flanks of the dome are vertical. In this case the area A 
is formed by two rectangles (see Fig. 5*) and there are seven parameters: the radii a 
(salt) and 6 (cap-rock) of two cylinders which together form the dome with its cap- 
rock, the depth # at which the salt begins (bottom of the cap-rock), the thickness h 
of the cap-rock and that H of the 
p salt, the density-contrasts k; and 
ke (salt). The expression (31) gives 
the characteristic function W(z), 


a 4?W(u) = e~“? | akoJ 1(au) 
(1 — e-“#) + bhyJi(bu)(e — 1) ]. 


? Choosing any seven numerical val- 
ues for u and computing for them 
the corresponding values of the 
product ~!u?W(u) with the aid of 
the rule (38), we have a system of 
seven equations with seven un- 
knowns 4, b, ki, ke, p, H, h. Solving 
it we have all the necessary infor- 
mations about the dome. This 
method can be applied only if the 
flanks of the dome are vertical. Let 
us now consider the case of a 
deeply buried salt dome, substi- 
tuting for it as an idealized form a 
truncated cone with the vertex at O*. This assumption is not acceptable for shallow 
domes. We need the moment function and we begin by giving its general expression 
which holds for all solids of revolution. A. Erdelyi and H. Kober have recently 
proved? an important theorem on the Hankel transform which they formulate using 
Tricomi’s form of this transform. In our notation (Hankel’s form) this theorem 
states that if W=H,.2.w, where s>—1, a>0O, and u'/*w(u) belongs to L.(0, ~) 
then 7,.(W) =H,|T,.(w)], the operator 7,. being defined by 








CAP ROC 


SALT 


|= 


$ 
Fic. 5. 


P(a)Traf(2) = 2a" f (98 — 2)ty44(y)dy. 
Applying it to our relation W(u) =H )Dg(r) with —1<s <0, s+2a=0, we have 


J ot - sremmmoyyay = +f seaneenan fe — eee 
z 0 


* In Fig. 5 the letters A; and A; denote the centers of the two end surfaces of the cap rock cylinder. 
2 Quarterly Journal of Math., 2, 217 (1940), Theorem 5, case b. 
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Interchanging the integrations in the second member (absolutely convergent double 


integral), using t=rv'/? and passing to the limit x =0, we get the desired general ex- 
pression for the moment function M(s) of Dg(r) for —2<s <0, 


M(s) = Sf. Dg(r)r'dS = 2 f° Deeper 


22+eT2(1 + 5/2) sin (— sx/2) f W(wu-rdu, (40) 
0 


The same result can be obtained integrating the relation r*+!Dg = Hor*t!W(u) un- 
der the sign of integration in u and applying the integral (33). Now the integral in 
the second member of (40) can be calculated, using (31) and the formula (3) p. 385, 
ch. 13.2 of “Bessel Functions” by G. N. Watson, 1922. Let w be the angle the radius 
vector of a point (p, ¢) makes with the z-axis; thus { = (p?+£*)"/? cos w, p=¢ tan w. Let 
P,, be the Legendre function of the first kind P,(cos w). Since F[(1+s)/2, —s/2; 1; 
sin? w] =P, and (2+s) sin? wF[(3+5)/2, —s/2; 2; sin? w] =2(P,—cos wP,4:), we have 


M(s) = (2+ s)e, ff k sec* wP,f*pdpdt = af k sec?+* w(P, — cos wP,4:)f**+*dt, (41) 
A c 


where (2+5)c, =2rT(1/2)T'(1+s/2)T(1—s/2). The formula (41) solves the problem 
for a deeply buried dome. The parameters in this case (see Fig. 6) are: the three 
lengths O*A ;=p;, j=1, 2, 3, the density contrasts ki =2foi, ke=2fo2 and the angle w, 


Oo 





| vm 























Fic. 6. 


the slope of the dome’s flanks. Instead of p; we use the ratios r= pi/p2, g=ps/p2 and 
write p,=>. Also, k=ke/k; this ratio being generally negative. The parameters r<1, 
q>1, k verify the equation 

4c*[1 — r? + (g? — 1)k] = 3[1 — r* + (qt — 1)k]p. (42) 


In our case w is constant in the second integral of (41) since only the integral 








74 E. G. KOGBETLIANTZ [Vol. III, No. 1 


along BA does not vanish. It is equal to (3+s)kip*+*o(s; r, g, k), where (s; 1, g, R) 
=1—r3+*+h(g3+*—1). Denoting (3+5)~'c, sec?+* w[P,(cos w) —cos wP,4:1(cos w)] by 
f(s; @), we have the explicit expression of the moment function M(s) 


M(s) = kip***f(s; w)@(s; 7, 9, ). (43) 


The f(s; w) is easy to tabulate for eleven particular values of s, namely for 4s =m, 
where the integer m verifies —7 Sm 33 since —2<s<1. If 4s =m the function P, re- 
duces to polynomials or to elliptic integrals K and £. In fact, letting cos w=x, we 
have 3P_7/4=3P3j4=Pijat+2xP_ijs, P_s a= Priya, P_s; a= Piya, SP 5/4= 6xP a— Piya, 
21P1j4=10xP1j4+(20x?—9)P_ij4, where m cos (w/4)P_1;4=2[cos w/2)]'4K and 
m(cos (w/2)]'/4*Pij4 = 2 cos w/4)[4(E — K) + (1+22)K] with the modulus 
A=tan (w/4). 

Also P_32=Pij2, 3P32=4xP12—P_1/2, where P12 are elliptic integrals too, but 
with modulus u=sin (w/2), namely tP_12=2K and rPi.2=4E—2K. Transforming 
the experimental gravity map into the maps of the quantities r*Dg(r) and integrating 
on these maps over r SR, we obtain the reduced moment function M*(s) which differs 
from M(s) only in the contribution of the infinite area r2R. This contribution is 
computed with the aid of (36), and the reduction factor p(s) in 


M(s) = M*(s)[1 + »(s)(p/R)'-*] (44) 


is defined by 4(1—s)f(s; w)(s; 7, g, k)v(s) =3f(0; w)d(1; 7, q, R). 

We consider equations of the general type Q(s, t; w, r, g, k) = N(s, t), where the 
function Q of four unknowns a, r, q, k is defined by Q=[f(s; w)o(s) |-*[ f(t; w)o(d) |* 
- [f(0; w)(0) |*-*, the o(s) denoting $(s; r, g, k). Thus the number N(s, ¢) is to be 
calculated by the rule N(s, t) =[M(s) ]-*[M(2)]*[M(0) |*-*. Its first value is obtained 
by using the reduced moment functions M*(s), M*(t), M@*(0). With the aid of these first 
values N* we have to find first approximations for our parameters w, r, g, k. It is 
sufficient to form four equations of the type Q=N, giving to the orders s and ¢ nu- 
merical values, for instance s=—7/4 and t=—5/4, —3/2 and —1, —5/4 and 
— 3/4, —1 and —1/2. Solving such a system with the experimental data N*(—7/4, 
—5/4), N*(—3/2, —1), N*(—5/4, —3/4) and N*(—1, —1/2), we get the first ap- 
proximate values for w, r, g, k. The depth p is then found using the ratio M*(s)/M*(0) 
for any value of s (or ¢) and in practice the average value from many such determina- 
tions will be taken as p. Using the set of first approximations in (44), we improve the 
first values of the second members Ms, ¢) in our system of four equations Q = N and 
solving the improved system, we have second approximations for w, r, g, k, p. This 
procedure of successive approximations is continued until the stabilization of se- 
quences. The depths p:=rp and p3;=qp are found, as well as ke=kki, since the value 
of ki can be obtained from the value of any M(s). We observe that the interpretation 
yields both density contrasts o; (cap-rock) and g (salt). The depth 2* of the center of 
gravity obtained with the aid of (42), compared with z* computed by (23*), gives a 
control. The control also is obtained with the aid of four other values of the order s 
which were not used for the interpretation. We choose the negative values of s in 
order to diminish the reduction factor in (44). It is interesting to add that the same 
method applies to the interpretation of G- and K-maps, their moment functions 
G(s), K(s) being expressed in terms of M(s). In fact G(s) =—(2+s)M(s—1) and 
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K(s) =c} M(s—1) with c/ T(1+s/2)T(1/2 —s/2) =21T(—s/2)T(5/2+s/2). For the in- 
terpretation of the magnetic maps we can use the same function M(s) since the 
moment functions A,(s) and Z(r) are given by ke(2+s)Z(s) = —2I sin y-sK(s) and 
koA ,(s) = —2] sin Y-(2+5)M(s). 

In the general case, when the vertex of the cone is on the vertical of the point O* 
at the distance (unknown) / above (J>0) or below (<0) the ground, the moment 
function has a very complicated expression and is difficult to tabulate. In its place 
we use the integrals D, of the type 


D,, = pe f f (L? + r?)— (+112) g(r) dS, 
P 


and the interpretation yields all the seven parameters hk, ke, pi, p2, ps, w, 1, the angle w 
being the slope of the dome’s flanks. Lack of space does not permit the development of 
this general case. 

Conclusion. The possibilities offered by the new method for quantitative inter- 
pretation of magnetic and gravitational anomalies we attempted to develop in this 
work seem to be very large. The mathematical tools, used in the proofs of the final 
interpretation rules, are not needed at all in practical applications. If the tables of 
functions used and the charts of auxiliary curves for the graphical solution of funda- 
mental equations are calculated once for all and plotted, all that remains to be done 
in a particular case is the plotting of some auxiliary maps derived from the experi- 
mental data and the evaluation of some areas with the aid of a planimeter. 

No use of average values in the interpretation is known to the author, with the 
exception of the work of K. Jung,* where only one integral, namely our M(0), is 
defined and its value 2fM is used together with remarkable values and distances. 


3K. Jung, Zeit. f. Geophysik 13, 45-67 (1937). 











76 


CANTILEVER BEAMS OF UNIFORM STRENGTH* 


BY 
I, OPATOWSKI 


Armour Research Foundation 


1. The object of the paper, its methods and results. The problem of shaping a 
beam from a given amount of material in such a manner as to obtain maximum 
strength requires that the maximum stress of each cross section be constant. In the 
case of bending, the classical treatment of this problem!*-4 is based on the theory of 
beams of constant cross section, the influence of shearing stresses and of the weight 
of the beam being neglected. A collection of solutions of this elementary problem, 
for rectangular and circular cross sections, is given in the Hiitte handbook for engi- 
neers.® If the strength of the material is relatively low, the weight W of the beam 
cannot be neglected. This occurs in certain concrete structures, such as reinforced 
concrete bridges, and was demonstrated by Gaede' in his treatment of a cantilever 
of rectangular cross section and constant width, the external load being a force F at 
the free end. 

In the present paper, we shall consider cantilevers of more general cross section 
but with the same type of loading, except in §6 where more general loading will be 
considered. Let us denote by x the distance from the free end, by A (x) the area of the 
cross section and by S(x) the section modulus (S= M/o, where M is the bending mo- 
ment and @ is the maximum stress). The bending moment M(x) =aS(x) at the dis- 
tance x from the free end is then given by 


Fxt+ vf (x — E)A(E)\dE = oS(x) (1.1) 


0 


where y denotes the density of the beam material. The total weight of the beam equals 
L 
vf A(t)dt = W (1.2) 
0 


where L is the length of the beam. Since ¢ is constant along the beam, differentiation 
of (1.1) with respect to x yields 


* Received March 20, 1944. This work has begun at the University of Minnesota and was completed 
at the Armour Research Foundation. Presented to the American Mathematical Society under different 
titles at the meetings of September 12-13, 1943 and February 25-26, 1944. The author is indebted to 
Professor G. E. Hay for many valuable improvements which were included in the text of this paper. 

1 F, Grashof, Theorie der Elasticitat und Festigkeit, R. Gaertner, Berlin, 2nd Edition, 1878, pp. 113- 
121. 

2 C. Bach, Elasticitat und Festigkeit, J. Springer, Berlin, 2nd Edition, 1894, pp. 85-88. 

3S. Timoshenko, Strength of Materials, part I, 2nd Edition, D. Van Nostrand Company, New York, 
1940, pp. 209-210. 

4 C. Guidi, Teoria dell’ elasticita e resistenza dei materiali, 11th Edition, Torino, 1925, pp. 135-142. 

5 Hiitte—Des Ingenieurs Taschenbuch, vol. 1, 25th Edition, Wilhelm Ernst, Berlin, 1925, pp. 626-629. 

6 Gaede, Balkentrager von gleichem Widerstande gegen Biegung, Bautechnik, 15, 120-122 (1937). 
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aS’(x), S(0) = 0, (1.1’) 


F+ vf A@at 
yA(x) = oS"'(x), S(0) = 0, oS'(0) = F, (1.1”) 


where the primes denote derivatives with respect to x. By use of (1.1’) we can write 


(1.2) in the form 
oS'(L) = F + W. (1. 2") 
We note that (1.1’) and (1.1’’) are forms of the well-known equations of equilibrium 
of a beam, Q=M’, g=M"’, where Q, g are respectively the shearing force and load 
per unit length. 

If the section modulus is assigned, A(x) is given by (1.1’’) and the problem is 
solved. In general however there are no criteria for the choice of the function S(x); 
instead, some geometric characteristics of the cross section are assigned. Problems 
of this type are treated in the present paper in a general manner. They involve an 
integral equation (cf. Blasius’). Its solution may involve almost any of the classical 
special functions. Some simple cases leading to hyperbolic, Bessel and elliptic func- 
tions are discussed. The possibility of using Legendre, hypergeometric, Lamé and 
some other functions is indicated. 

2. The type of beam. Throughout this paper we shall limit ourselves to cantilevers 
satisfying the following conditions: the line of centroids is a horizontal straight line 
(x-axis); each cross section has a vertical axis of symmetry (V-axis). In the plane of 
the cross section we choose a system of orthogonal Cartesian coordinates (U, V) with 
origin at the centroid C. In the vertical plane through the x-axis, we choose a system 
of Cartesian coordinates (x, y) with origin at the free end and y-axis directed down- 
ward. We assume that the curves bounding the cross sections are representable by 


the equations 
U = u(x)u,(2), V = v(x)0,(2), (2.1) 
t being a parameter. The functions ,(¢), v:(¢) determine the shape of the cross sec- 
tion, whereas the functions u(x) and v(x) represent the change of the cross section 
along the axis of the beam. Any two cross sections are obtainable from each other by 
a transformation of dilatation* which depends on the position of the cross sections. 
We will choose u(t) and »;(t) in such a manner that u(x) and v(x) be 20. 

3. General equations. It is easily seen that S=I/V», where J is the moment of 
inertia of the cross section about the U-axis, and V, is the maximum value of V. 
Thus, if @ is the area enclosed by the curve U=1,(t), V =2;(t) and 8 the corresponding 


section modulus, we have 
S = Bu(x) [v(x) ]?, A = au(x)v(x). (3.1) 
If we set a=ay/(o8), the substitution of (3.1) into (1.1), (1.1’’), (1.2), (1.2’) gives 


Fat ay f(x — du(do(@de = o8u(2)[o(2) (3.2) 
0 
7 H. Blasius, Trager kleinster Durchbiegung und Stabe grisster Knickfestigkeit bei gegebenem Materialver- 


brauch, Zeit. f. Math. u. Phys., 62, 182-197 (1914). 
8 The writer is indebted to Prof. H. Busemann for this geometric terminology. 
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(uv®)”” = auv, (uv?) .-0 = 0, oB(uv?) 0 = F, (3.7) 
L 
« f u(x)o(x)dx = W/y, (3.3); oB(uv’) 1 =F + W. (3. 3’) 
0 


If v(x) is known, (3.2) is a Volterra integral equation in u(x) with the kernel 
(x —£)/[v(x) ?. This kernel is a continuous function, within the interval of integration, 
if v(0) +0, because we assume v(x) continuous and by its physical meaning it must 
be +0 for x>0. Therefore, according to the general theory of integral equations,? if 
v(0) #0, Eq. (3.2) has one and only one solution u(x) if F¥0 and only a meaningless 
solution «=0 if F=0. In other words, a cantilever of uniform strength under the ac- 
tion of its own weight alone must be such that v(0) =0. 

4. Particular types of cantilevers of uniform strength. These are obtained by as- 
suming particular forms for u(x) or v(x). 

I). v is constant. The cross sections have constant height. If F+0 the integral of 


(3.2’) is 
u = Fr(ayv)— sinh (rx) (4.1) 


where 7 = (a/v)!/?. Substitution from this equation in (3.3’) gives the following condi- 


tion for W: 
cosh (rL) = 1+ (W/F). (4.2) 


If F=0 no solution exists, which is in accordance with the general statement of §3, 


because here v(0) 0. 
II). v(x) is a linear function of x. The cross sections have linearly varying height. 


We may restrict ourselves to the case 
v(x) =ctx (4.3) 
since if x had a coefficient different from +1, the coefficient could be factored out and 


included in the function 2;(¢). Also, since we agreed to take v(x) 20 (cf. §2) and x=0 
represents a point of the beam, c must be 20. Since dv = +dx, the solution of (3.2’), 


(3.3) is? 
u = v-¥/2Z,(2ia'%y'/2), (4.4) 
where Z, is a cylindrical function of order 1 which must satisfy the conditions 
Z,(2ia'/*¢1/2) = 0, + oBa'!*iZo(2ia/*c1/?) = F, (4.5) 
+ oBa'!iZo[2ia/*(c + L)'/2?] =F + W. (4.5’) 


The second equation in (4.5) and Eq. (4.5’) are obtained by use of the formula’ 
Zi (x) =Zo(x) —x1Z;(x). We put 


1/2 1/2 1/2 


Z,(2ia o ) = AiJ;(2ia » ) + BH; (2ia’ 0 


{2 


), (4.6) 


where J; and H!” are the Bessel and the Hankel*functions of the first kind and first 
order.!° Equations (4.5), (4.5’) then give the following conditions for A, B and c: 


® See for instance R. Courant-D. Hilbert, Methoden der mathematischen Physik, vol. 1, 2nd Edition, 
Julius Springer, Berlin, 1931, pp. 119, 120, 133; or Riemann-Weber, Die Differential- und Integral- 
gleichungen der Mechanik und Physik, vol. 1, 7th Edition, F. Vieweg, Braunschweig, 1925, pp. 426-428 
or E. Hellinger-O. Toeplitz, Enc. d. math. Wiss., vol. II. 3, pp. 1459, 1460. 

10 E, Jahnke-F. Emde, Funktionentafeln mit Formeln und Kurven, 3rd Edition, B. G. Teubner, 


Leipzig, 1938, pp. 128, 144-147, 224, 237, 242. 
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AiJ,(2ia"c!”) + BH (2%a'*c'”) = 0, (4.7’) 
— AJ((2ia’"c'”) + BiH, (2ia'c “ = +F/(oBa’), (4.7") 


+(F+W)/(oBa’). (4.8) 


We discuss first the case c=0, i.e., v(x) =x (the lower sign in (4.3) has no meaning 
here, since v must be 20). The free end of the cantilever is represented by v=x=0. 
Therefore, since H{”(0) = «, the constant B must be zero. Since J,(0) =0 and Jo(0) 
=1, Eqs. (4.7’), (4.7”) require that A = — F/(aBa"/?), and Eq. (4.8) gives 


Jo(2ia'/2L'/2) = 1+ (W/F). 


If L and W are not related by this equation, the constant c must be distinct from 
zero. The determinant of the coefficients of (4.7’), (4.7”), considered as equations 
in A and B is, by a known relation of Bessel functions,'° 


— AJo[2i(ac + aL)”] + BiH [2i(ac + aL)” 


(1) (1) —1 —1/2 —1/2 
Jo(z)Hy (2) — Ho (2)Ji(s) = —e ac , (4.9) 
where z= 2ia‘/*c!/*. Since this cannot be zero it is seen that there are no solutions if 
F=0, which agrees with the general result of §3. If F#0, the solutions of (4.7’), 
(4.7”) are 
1/2 u2 we 


A= +-2(08) ¢ FH, (Qia c’), B= (08) c FiJ\(Qia c). (4.10) 
Substitution from these into (4.8) gives a relation for W,™ 
‘c [Ji(z)Ho (&) — Hy (2)Jo(t)] = 1+ (W/F), (4.11) 


where 2=2i(ac)'/?, €=2t(actaL)"/?. 

III). u(x) is proportional to [v(x) |". This includes a circular cross section (n=1), 
a rectangular cross section of constant width (m=0), a rectangular cross section with 
the height proportional to the width (m=1), an elliptic cross section with axes pro- 
portional to each other (n=1). By a suitable choice of u(t) and 2;(¢) we may reduce 
the problem to the case u=v". The first two equations in (3.2’) then give 


+= i) v™t1[C2 + 2a(n + 2)-(2n + 3)—!02*+8 |-1/2dp, (4.12) 
0 
where C is a constant. The last equation of (3.2’) and Eq. (3.3’) give 
oBC = F/(n + 2), = [F’' + (n+ 2)(2n + 3) 2yoabv, 1) —F. (4.13) 
If m= —1 the cross sections have a constant area; this case gives elementary expres- 


sions for u(x) and v(x) but the width at the free end is infinite. When »=1, the in- 

1 By (4.3) the possible range of L is (0, + ©) or (0, c). The left hand side of (4.11) is within this range 
a function of ZL which increases continuously from 1 to +, as may be shown by the theory of Bessel 
functions and by (4.7’), (4.9), (4.10). Therefore for any ouslanall values of F, W, a, c there exists, in each 
case (4.3), one and ealy one value of L satisfying (4.11). Problems of this type in which the length L is 
not assigned but a given amount of material is to be distributed into a cantilever of uniform strength 
under the action of F and W or W alone have occurred in some biological fields (cf. C. Holtermann, 
Schwendener’s Vorlesungen iiber mechanische Probleme der Botanik, Leipzig, 1909, pp. 18, 19 and O. Fischer, 


Enc. d. math. Wiss., 1V.8, p. 119). 
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tegral in (4.12) is hyperelliptic; when n =0, it is elliptic. If F=0, Eqs. (4.12) and (4.13) 

give 

vy = a[2(2n + 3)(m + 2) }-'x?. (4.14) 

IV. More general cases. Instead of u and v we introduce new variables w= uv’, 

Tt =1/v; is directly proportional to the section modulus, and 7 inversely proportional 
to the radius of gyration of the cross section. From (3.2) we obtain 


Pid f (x — £)r()w(8)d = oBwo(2). (4.15) 


This is a Volterra integral equation in w(x). From it, or directly from (3.2’), (3.3’) 


we have 
w’’(x) = ar(x)w(x), w(0) = 0, ocBw'(0) =F, cBw'(L) =F + W. (4.16) 


Since most of the so called special functions satisfy linear differential equations of 
second order, the first equation of (4.16) suggests the possibility of using such func- 
tions. The following are some results which may be easily checked. The constants 
P, 7, Ss, @ must satisfy the last three equations in (4.16). 

IVa). r(x) =p—ge**, w(x) =Z,,(ne*), where Z,,=a cylindrical function (Bessel, 
Hankel, etc.), m?=ap, n?=aq. } 

IVb). r(x) =p—q(cosh x)-?, w(x)=K” (tanh x), where K®=an associate Le- 
gendre function (P*’, oO”), m?=ap, n(n+1) = aq. 

IVc). r(x) =p—gq cos x, w(x) =a function of an elliptic cylinder.” 

IVd). r(x) =(p—qx+x?)/(4ax?), w(x) =a confluent hypergeometric function. !* 

IVe). r(x) =(p—qx*)/x?, w(x) =x"/?Z,,(nx*/*), where Z,,=a cylindrical function, !® 
m*s? =1+4ap, n*s*=4ag. If p=0, in order that v be finite s must be <2. 

If the function v(x) =1/r(x) is assigned by means of any one of previous expres- 
sions for 7, the function “=wr? is determined by the corresponding expression for 
w(x). In the case of a rectangular cross section, v(x) represents the height and u(x) the 
width. 

5. The deflection curve. The curvature of the geometric axis of a beam of con- 
stant strength in bending is!-? 1/r =h/v(x) where h=a/(Ev,,), E being the modulus of 
elasticity and v,, the value of v,(¢) at the point of maximum stress (cf. §2). We note 
that this equation is a form of the well-known relation ¢ = Ey/r. For small deflections 
the usual approximation is 1/r =d?y/dx?. Thus 


z L 
y(x) = -f g(x)dx, g(x) = nf [v(x) |-'dax, (5.1) 


since 
‘ y(0)=0, (dy/dx) 1 = 0. (5.2) 


A simple formula for the deflection at the free end is obtained through integration 


of (5.1) by parts. Setting —y(L) = Y, we have 


Y=h {= {Ios} eee + h f° [o(2)-*2dz = bf [o(2) Fade (5.3) 
z / t=0 0 0 


2 E, T. Whittaker-G. N. Watson, Modern analysis, 4th Edition, Cambridge University Press, 1935, 
pp. 337, 405. 
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It is seen from (5.3) that, if v(x) tends to zero as kx® with n=2 and k is constant, the 
deflection Y is infinite. This would occur, for instance, in the case corresponding to 
Eq. (4.14). Such a physically impossible conclusion may be explained by the fact that 
a large value of » implies a rapid variation of v(x), i.e., a rapid change of the cross 
section, whereas the theory which was used is based on bending of beams of constant 
cross section."* More important still, the theory used in this paper neglects the shear- 
ing stresses in comparison with the bending stresses. Such a procedure is not permis- 
sible in the vicinity of the free end, and consequently it is understandable that the 
theoretical results for this part of the beam differ widely from reality. 

6. More general loads. If M(x) is the moment of the external load acting on the 
cantilever, we have instead of Eqs. (1.1’’), (1.2’) 


M(x) + yA(x) = oS"(x), oS(0) = M(0), —o S’(0) = M"(0), (6.1) 
oS'(L) = M'"(L) + W. (6.2) 
For example, if the beam is acted upon by F and also by a load distributed uniformly 


along the axis of the beam of intensity T, we have M(x) = Fx+4Tx?. If v=const., 
we obtain by (6.1), (3.1), and (6.2), 


u = [Fr sinh (rx) + T cosh (rx) — T]/(ey») (6.3) 
cosh (rL) + T(Fr)“ sinh (rL) = 1+ (W + TL)/F, (6.4) 


where r=(a/v)"/*. If T=0, Eqs. (6.3), (6.4) reduce to (4.1), (4.2). Eqs. (6.3), (6.4) 
may be easily generalized to the case M(x) =)>a,x". 

7. Numerical examples. We consider a rectangular cross section of width 1 and 
height H. Then (cf. §3) a=H, v,=H/2, B=H?/6, a=6y(cH)—'. Let L=10 ft., 
F=9000 Ibs., ¢ = 75000 Ibs./sq. ft., y=150 Ibs./cu. ft., E=45 X10" Ibs./sq. ft. These 
values correspond to a certain type of concrete. 

I). Cantilever of constant height (§4, Case I). We put H=1 and assume the height 
vH =v=1.9 ft. From (4.2) we obtain the weight W = 3000 Ibs. We put R= [6y/(ov) ]"/*. 
Then R=0.0795. From (4.1) we obtain the width 


u(x) = FR(vy) sinh (Rx) = 2.51 sinh (0.0795x). 


At the fixed end we then have u=u(10) ~2.21 ft. Equation (5.3) gives for the deflec- 
tion Y=oL?/Ev~0.1 in. 

II). Cantilever with a linearly varying height (§4, Case II). Let the height at the 
fixed end be 2 ft. and at the free end 1/4 ft. In (4.3) we take v(x)=c+-x. Since 
Hce=1/4 ft., H(c+10) =2 ft., we get H=7/40, c=10/7 ft. Eqs. (4.10) give A = —69.1, 
B=29.0, and from Eqs. (4.4), (4.6) we obtain 


u=v  [— 69.14J,(it) + 29.0H, (it)], where £=0.8(30/7). (7.1) 
At the fixed end x=10, and we thus obtain u=2.,2 ft. At the free end, v=c and by 
(4.4), (4.5) we have the general result u=0. From (4.5’) or (4.11) we get the weight 
W ~4800 Ibs. From (5.3) the deflection is 
Y = 20(EH)™|L — clog. (1 + Lc)] = 0.2 in. 


18 A method which takes into account the variability of the cross section was worked out by J. Résal, 
Résistance des matériaux, Paris, 1898, pp. 393-405 for rectangular and double T cross sections of con- 


stant width. 
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THE INVERSE OF A STIFFNESS MATRIX* 
By K. E. BISSHOPP (Armour Research Foundation) 


It is well known that torsional vibration problems often require the computation 
of latent roots of matrices. Now the usual methods! give these roots in descending 
order of magnitude while in torsional vibration problems we require the smallest 
root of the stiffness matrix and then, perhaps, some of the remaining roots in ascend- 
ing order of magnitude. It is therefore necessary to find the inverse u—' of the given 
stiffness matrix; u~ is called the flexibility matrix. If its roots are in descending order 
pi, p2, ++ +, Px then the required roots of the original stiffness matrix in ascending 
order are 1/1, 1/po, +--+, 1/px. 

In general it is very difficult to invert-a given matrix. The purpose of this note is 
to show that a special type of stiffness matrix which occurs frequently can be inverted 
with a small amount of work. Let us consider, for purposes of illustration, +1 discs 
with moments of inertia Jo, - - - , J, connected by massless elastic shafts of circular** 
cross section. Let c; be the coefficient of stiffness of the shaft between the zth and 
(t+1)th discs. Then 

c; = GJ,/li, (1) 


where /; is the length of the shaft in question and GJ, is a numerical factor depending 
on the material and the polar sectional moment of inertia. The application of La- 
grange’s equations? to the functions representing the kinetic and potential energies 
of the system respectively yields a system of linear differential equations. Therefore 
the second order time derivatives can be replaced by — pj, so that in the absence of 
damping we obtain the following system of algebraic equations in matrix form, 


(p> — u)6 = 0, (2) 


where 8 is the column matrix of the normal mode appropriate to any value of p? for 
which Eq. (2) is satisfied. Since the system is capable of a rigid body rotation, p?=0 
is a solution and the degree of system (2) can be reduced by unity with the aid of the 


substitution 
To + I 10; + tone: + TOn = 0, 


* From a paper, The use of matrices and normal coordinates in the solution of torsional vibration prob- 
lems, read at the spring meeting of the Wisconsin section of the Math. Assoc. of America, at Milwaukee, 
Wis., May 13, 1944. Manuscript received August 18, 1944. 

1 R. A. Frazer, W. J. Duncan, and A. R. Collar, Elementary matrices and some applications to dynamics 
and differential equations, Cambridge University Press, Cambridge, 1st Ed., 1938, Chapters I-V. 

** For other than circular cross sections, suitable factors can be obtained from the appropriate torsion 
functions. 

2S. Timoshenko, Vibration problems in engineering, 1st. Ed., D. Van Nostrand Co., New York, 1928, 


p. 139. 
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which may be looked upon as the orthogonality condition*® for the degenerate fre- 
quency p=0. The stiffness matrix u for the reduced system in which the zero root is 


absent becomes 














Co Co Cy Ci Cole Cols Coln a 
— a ints jin sonepeniinne ese 
I I; q; I, Tol Tol, Tol, 
Cj Ci C2 C2 
a pee. = ae 68 tent ose 0 
I, Ie I. I. 
Ce Cc Cc : 3 
4 ss a). (3) 
TI; I; Is 
0 0 0 a 
b eS 








The inverse matrix u~!=(a;;) may be stated in a convenient form for numerical com- 


putation as follows. 
If A =>" J,, then the jth element in the first row becomes 


ans = [GE DI (Iolo — L tetvss). (4) 


When i Sj the ith element in the jth column is 


i-1 
aij = 1;+ (1 ;/GI«) >, les (S) 
s=1 
and when 7>7, 
aij = Q;;- (6) 


As an application, let us consider the case when Jo = 181 . 306 Ib.in.sec.*, J; =J2=Ts 
= 1,328.61 Ib.in.sec.2, J4=21,557.3 lb.in.sec.?, 19 =30 in., 1; =13 = 34 in., 13 =62.2 in., 











TABLE I 
(1) (2) | (3) (4) | (5) 

. Be lL. I, X10°/GJ, An | Ido oP nA ngs 
0 181.306 30. | 25,724.4 

1 1,328.61 a 0.076051 | 25,543.1 | 5,439.2 
2 1,328.61 34. 0.076051 | 24,214.5 | —817,850. 
3 1,328.61 62.2 0.076051 | 22,885.9 | —1,595,970. 
4 21,557.3 1.23396 21,557.3 | —2,936,840. 





GJ, =17,470 X 10° Ib.in.? All necessary computations are contained in Table I. From 
the second line of this table (n=1), we obtain immediately 


10®a1, = Iolol1/GJ,Ao = [col. (5) X col. (3)]/Ao. 


3T. von Kérm4n and M. A. Biot, Mathematical methods in engineering, McGraw-Hill Book Co., 
New York, 1940, chapter V, pp. 162-215. 
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Equation (4) shows that the remainder of the elements of the first row of u~' can be 
calculated in order from each succeeding line of the table by performing the operations 
indicated on columns (3) and (5) respectively. The other elements in the upper right- 
hand corner of u~! then are computed from Eq. (5) which gives for instance when 
4=3 and j=4, 


10%a3, = 10%ay4 + Ty X 10°/GJ, >, Il, = — 140.875 + 1.23396 X (34 + 34) 


n=1 


= — 56.966. 


The inverse matrix now can be completed quickly by filling in the lower left-hand 
corner according to Eq. (6), so that 


0.0160803 — 2.4179 — 4.7183 — 140.875 
- 0.0160803 0.1678 — 2.1326 — 98.920 
oie 10-°*X] | 

0.0160803 0.1678 0.4532 — 56.966 

0.0160803 0.1678 0.4532 19.786 


ON THE PROBLEM OF HEAT CONDUCTION IN A 
SEMI-INFINITE RADIATING WIRE* 


By ARNOLD N. LOWAN (Math. Tables Project, Nat. Bureau of Standards) 


R. V. Churchill! derives the solution of the problem of heat conduction in a semi- 
infinite radiating wire when the initial temperature is zero, and the boundary tem- 
perature is a constant. It is the object of this paper to derive the general solution 
corresponding to an arbitrary initial temperature distribution when the boundary 
temperature is a prescribed function of time. 

Let k, c, p, s, A, h and a=k/pc denote the thermal conductivity, specific heat, 
density, perimeter, cross-sectional area, coefficient of heat transfer, and thermal diffu- 
sivity of the wire, respectively. Further, let a=hs/cpA and b=aT>2, where T:2 is the 
temperature of the medium. If the wire is sufficiently thin so that the temperature 
may be assumed to be constant over the entire cross section, the problem becomes 
one-dimensional and the temperature T(x, t) must satisfy the following differential 
equation, initial and boundary conditions: 





0 0? 
(< - a +a)r(x9 =3 (x > 0,¢#> 0), (1) 
yom T(x, t) = f(x), (2); T(0, t) = ¢(?). (3) 


It is easily verified that the expression 
T(x, t) = e~*'u(x, t) + 0(2, é) (4) 


* Received July 17, 1944. 
1R. V. Churchill, Modern operational methods in engineering, McGraw-Hill Book Company, New 


York, 1944, p. 119. 
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satisfies Eqs. (1), (2) and (3), provided the functions u(x, ¢) and v(x, ¢) satisfy the 
following differential equations, initial and boundary conditions: 





Ou O24 
—-—-a =x ( (x > 0,#> 0), (5) 
ot Ox? 
lim u(x, t) = f(x), (6); u(0, t) = e*g(t), (7) 
t-0 
Ov 02y 
ae av+b, (8); lim v(x, 2) = 0, (9); (0, 2) = 0. (10) 
€ Cos” 1-0 


From Eggs. (5), (6) and (7), it is clear that u(x, t) is the temperature in a semi- 
infinite solid initially at the temperature f(x) and with its bounding plane x =0 kept 
at the temperature e*‘g(¢). Using the expression of u(x, t) given by H. S. Carslaw,? 


we obtain 


—at 





e—*'y(x, 1) = — 


f s@feevtraa =i eW(2t)*/(4a0) | de 
0 


2V rat 
x ' ee 
+—— f§ oe —- ger ye. (11) 
2V ma 7 0 


With the aid of the identity 
/ 


at 





e7 (2-8) */hat 


J e~*8"t cos B(x — £)dB = 
0 2 
the first term of (11) may be written in the alternative form 
2 co) co) 
— cof f e~ a8" tf(£) sin Bx sin BitdBdé. 
T 0 0 


Accordingly, an alternative form of (11) is 








Jeet 2 -) 
e~*'4(z, i) = f f e~ 28" tf(£) sin Bx sin BédBdt 
T 0 0 
x t 
+ — J e-*9(t — nem? 40%y-8/2dn, (11’) 
2V ra Y 0 


We proceed to the solution of the system (8), (9), (10). The Laplace transform 
v*(x, p) of the function v(x, ¢) must satisfy the equations 


02y* b 
pce Sto wee o*(0, p) = 0. (13) 


Ox? a pa 





The solution of the system (12) and (13) is 


b are 
o*(x, p) = ——— (1 — FST) = o*(p) — 0° (p) (9) (14) 
- Tipe p(p + a) 
2H. S. Carslaw, Mathematical theory of the conduction of heat in solids, Macmillan and Co., London, 
1921, §81, p. 172. 
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whence, by a well-known theorem, 





t 
v(x, 4) = a(t) -f a(t — n)w(n)dn. (15) 
From 
%0 b b/1 1 
sds iin F cB gt 
0 p(p + a) a\p prta 
it follows that 

b 

o(f) = —(1 — &*). (16) 
a 


If in the known identity‘ 


- x 
f e 4/- endl ty-8/2g¢ = e-2VNP 
0 T 


we put \=x?/4a and replace p by p-+a, we obtain 


. 20 
x ° —_——— 
= A 2 ae : 
a -f e pte ate 2 4aty 3/2q} = eV (9 +a)/a : 
0 


2\/ Ta 
whence 
x ‘ 
ws, = —— e~ate—z"/daty—3/2 (17) 
2\/ Ta 


In view of (16) and (17), (15) becomes 
v(x, ) -~— -€ ne ;1 tah " lke St is dany—3/2dy 4+ — (1 on e~**), (18) 
2av Ta 0 a 


Making use of the identity® 


) VT 
f e— (4482) X27) — —— cosh 2ab 
. 


2a 


JT , b Jr b 
+ — e*Ert(— — ac) — — etErf(— + ac 
4a c 4a c 


and some elementary transformations, we may write (18) in the alternative form 


( b fi 5 ( x Jt 
v(x, t)-= — — ¢ Er Saat 
a \ 2Vat/} 


b f a =— pre: x 
— 2 cosh x — a e mwa eErf \/at a ae ) 
2a x 2V at 


; oe x 
— ere “Ext ( Va + rer , (18’) 
2V al 


* See, for example, J. R. Carson, Electric circuit theory and operational calculus, McGraw-Hill Book 
Co., New York, 1926, p. 41. 

4 J. R. Carson, loc. cit., p. 39. 

5 This is a slight generalization of Churchill’s formula (4), page 120. 
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where 


2 z 
Erf(x) = — f e*'dp. 
Jaw o 


It should be noted that in (18’) the function v(x, ¢) is expressed in terms of tabulated 


functions. 
The final solution of our problem is given by (4) in conjunction with (11) and (18) 


or (11’) and (18’). 


THE SPHERICAL GYROCOMPASS* 
By WALTER KOHN (University of Toronto) 


In the existing literature on gyroscopes! the theory of the gyrocompass is de- 
veloped for the case of a rotor whose ellipsoid of inertia is an ellipsoid of revolution. 
The mathematics of this treatment is somewhat involved and, in deducing the differ- 
ential equations of motion, approximations based on the smallness of the earth's 
angular velocity are made. In the present communication we shall treat a gyrocom- 
pass the rotor of which has a spherical ellipsoid of inertia. The motion of such a 
gyrocompass is, of course, covered by the more 
general theory usually given, but owing to the 
symmetry of the sphere this case allows a con- 
siderably simpler, separate treatment in which, 
moreover, no approximations are necessary. At 
the same time the essential features of gryo- 


® 


scopic motion are preserved. e 

The following system will serve as a simple 
model of a spherical gyrocompass. The rotor is a 
rigid homogeneous sphere rotating freely about 
a light axle which passes through its centre. The 
ends of this axle can slide in a smooth horizontal 
ring which is concentric with the rotor and 
rigidly attached to the earth. When the rotor is 
set in rapid revolution about its axle the latter 
executes oscillations about the meridian which 








will now be examined. Fic. 1. 


In the figure the right-handed unit triad ,i, j,k, 
which is fixed relative to the earth is defined as follows: O is the center of the rotor; 
k lies in the direction of the upward vertical; i lies along the meridian and points 
north; j, pointing west, completes the triad. The unit vector, a, lies along the axle 
of the gyrocompass and the unit vector, e (in the i, k plane), is parallel to the earths’ 
axis; thus the angle A, between i and e, is the latitude of the observer. 

* Received July 10, 1944. 

1 Cf. T. Levi Civita and U. Amaldi, Lezioni di meccanica razionale, vol. 2, Zanichelli, Bologna, 1927, 
pp. 191-195; or J. L. Synge and B. A. Griffith, Principles of mechanics, McGraw-Hill Book Co., New 


York, 1942, pp. 430-433. 
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It is clear that the couple exerted by the ring on the rotor must be of the form 
G = G(a Xk). 


Further, if A is the moment of inertia of the rotor, w its angular velocity and h its 
angular momentum, we have the relation 


G =h = Ao. 


Consequently 
o-a = 0, (1); o-k = 0. (2) 


Since the gyrocompass has only two degrees of freedom, equations (1) and (2), to- 
gether with initial conditions, completely determine its motion. 

We observe that @ is made up of three parts: the spin of the sphere about its axle; 
the rotation of the axle relative to the frame i, j, k; and finally the absolute rotation 
of i, j, k or of the earth to which it is attached. Therefore we may write :w =sa+6k-+ Ne, 
where s is the spin of the rotor, 8 the angle between the meridian i and the axle a, 
and Q the angular velocity of the earth. Differentiation of this relation gives 
© = §a+sa+6k+6k, and since the angular velocity of a is 6k+Qe and that of k 
is Ne, this equation becomes 


® = a+ s(6k + Qe) X a+ bk + OQ(e Xk). (3) 


Substituting (3) into (1) and (2), we immediately arrive at 


$ — 62 cos Xd sin 6 = 0, (4); sQ cos sin6@+ 6 =0 (5) 


as the required equations of motion. 
The solution of these equations may be obtained in the usual way. From (4) it 


follows that s=s9+ cos A(cos 6)—cos 8), where so and 4 are the initial values of s 
and @. Inserting this value of s into (5), we obtain a differential equation for @ alone, 
which is of the classical type 6=f(@); this can be solved in terms of hyperelliptic func- 
tions. If the initial spin sp is great we may replace s in (5) by So to obtain the well 
known result : the motion of the axle a is identical with the motion of a simple pendu- 
lum, the position of equilibrium being in the direction of the meridian. 

An interesting property of the spherical gyrocompass can be deduced directly 
from Eqs. (4) and (5) which are exact. For, if we multiply by s and 6 respectively and 
add, we obtain ss=66=0, which on integration becomes s?-++6?=constant.? This 
shows that a spherical gyrocompass has an angular velocity of strictly constant 


magnitude relative to the earth. 
The author is indebted to Prof. A. Weinstein for his advice and criticism. 


2 Levi Civita and Amaldi, I.c., derive a corresponding result, namely 
346? + 4Cs? = const. 


for the general gyrocompass. (A and C are the transverse and axial moments of inertia respectively.) 
They then explain it by energy considerations. In fact, however, this equation and therefore also their 
reasoning is not strictly accurate. The exact form is 

3462 + 4Cs? + 3(C — A)Q? cos? d sin? @ = const. 
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